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1 Introduction 



The study of the effects of noise and fluctuations is a well established subject in sev- 
eral different disciplines ranging from pure mathematics (stochastic processes) to physics 
(fluctuations) and electrical engineering (noise and radiophysics). In traditional statisti- 
cal physics, fluctuations are of thermal origin giving rise to small departures from a mean 
value. They tend to zero as one approaches the thermodynamic limit in which different 
statistical descriptions (different ensembles) become equivalent. Likewise, in more applied 
contexts fluctuations or noise are usually regarded as small corrections to a deterministic 
(noise free) behavior that degrades a signal-to-noise ratio or can cause transmission errors. 
In such framework fluctuations are a correction that can be usually dealt with through 
some sort of linearization of dynamics around a mean or noise free dynamics. A different 
point of view about fluctuations emerges, for example, in the study of critical phenomena 
in the 1970's. The statistical physics description of these phenomena requires a formulation 
appropriate for a system dominated by fluctuations and nonlinearities. A linear theory only 
identifies the existence of a critical point by a divergence of fluctuations. 

The renewed interest and activity of the last 15-20 years on stochastic phenomena and 
their applications is precisely in the context of the study of nonlinear dynamics and insta- 
bilities in systems away from equilibrium. This activity has led to some new conceptual 
developments, applications, and new or rediscovered methodology |2|, |3|, |^, |5|, |6[. Among 
these we would like to emphasize here two very general aspects. One is the concept that 
noise need not only be a nuisance that spoils the "true" and desired behavior of the sys- 
tem, but rather noise might make possible new states and forms of behavior which do not 
appear in a noise- free limit. These situations might occur when there are mechanisms of 
noise amplification and/or when noise interacts with nonlinearities or other driving forces 
of a system. Phenomena like noise sustained spatial structures, noise-induced transitions 
or stochastic resonance go under this category. A second concept we wish to emphasize 
is that the physical relevant behavior is not necessarily associated with some ensemble 
average, but rather with typical stochastic trajectories. It is certainly a trivial mathemat- 
ical statement that a mean does not always give a typical characterization of a process, 
but there is a certain tradition in physics (borrowed from equilibrium statistical physics) 
of focusing on averaged values. A physical intuition or understanding of novel stochastic 
driven phenomena is in fact gained by the consideration of the individual realizations of 
the stochastic process. This has important methodological consequences: one needs tools 
to follow trajectories instead of just looking at averages and probability distributions. 

In these lectures we follow an incomplete random walk on the phase space of some of 
the current concepts, developments and applications of stochastic processes from the point 
of view of the practitioner physicist and emphasizing examples of the two key ideas out- 
lined above. Section 2 is rather tutorial while the other ones give more a summary and 
guide to different topics. Some parts of Sect. 2 are rather elementary and can be skipped 
by anyone with a little experience in stochastic processes. But this section also contains 
a rather detailed presentation of numerical methods for the simulation of the individual 
realizations of a stochastic process. It is here important to warn against the naive idea 
that including noise in the simulation of a nonlinear dynamical problem is just to add some 
random numbers in any reasonable way. This is particularly important when we want to 
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learn on the physics of the problem following the trajectories. Section 3 deals with one of 
the important cases of noise amplification, namely the transient decay from unstable states. 
Key concepts like trajectory dynamics, passage time statistics and mapping of linear into 
nonlinear stochastic properties are discussed. As an important application of these ideas 
and methods we discuss properties of laser switch-on viewed as a process of noise amplifica- 
tion. Section 4 is devoted to the analysis of the long-time, stationary, properties of physical 
systems in the presence of noise. The classification of potential/non-potential systems is 
succinctly reviewed and the conditions for a system to be potential are precisely stated. 
We study the general form of the probability distribution function when noise is present. 
We end with an example borrowed from fluid dynamics, the Kiippers-Lortz instability, in 
which we illustrate the stabilization by noise of a periodic trajectory in a system in which 
the deterministic dynamics has a contibution which can not be described as relaxation 
dynamics in a potential. In Section 5 we consider spatially extended systems, described 
either by ordinary differential equations for the amplitudes of a few spatial modes or by 
stochastic partial differential equations. Through some specific examples we discuss sym- 
metry breaking by pattern formation and symmetry restoring by noise, the issue of pattern 
selection in the presence of noise and noise sustained structures in convectively unstable 
situations. Section 6 reviews the concept of noise-induced transition, distinguishing it from 
that of noise-induced phase-transition. The difference being, mainly, that noise-induced 
transitions do not break ergodicity, as we understand a phase transition in the statistical- 
mechanics sense. We present specific examples of systems displaying one or the other and 
show that, in general, they can no coexist in the same system. 

We should finally make clear that our choice of subjects included in these lectures is 
rather arbitrary and sometimes dictated by personal contributions. Among many other 
relevant topics of actual interest which we do not discuss here we could mention, from the 
methodological point of view, path integral formulations of stochastic processes ^ and 
from the conceptual point of view stochastic resonance and directed fluxes supported by 
noise Q 



2 Stochastic Processes 
2.1 Basic Concepts 

In this flrst subsection we want to give a quick review of what a stochastic process is from 
the physical point of view. We will not be too rigorous on the mathematical side. The 
name "stochastic process" is usually associated with a trajectory which is random enough 
to demand a probabilistic description. Of course, the paradigmatic example being that of 
Brownian motion |^ 10, 11, 12, 13 1. The botanist Robert Brown discovered in 1827 that 



particles of pollen in suspension execute random movements which he interpreted initially 
as some sort of life. It is not so well known that L. Boltzmann knew as early as 1896 the 
reason for this erratic movement when he wrote "... very small particles in a gas execute 
motions which result from the fact that the pressure on the surface of the particles may 



fluctuate" [14|. However, it was A. Einstein in 1905 who successfully presented the correct 
description of the erratic movement of the Brownian particles ||2|. Instead of focusing on 
the trajectory of a single particle, Einstein derived a probabilistic description valid for an 
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ensemble of Brownian particles. In his description, no attempt is made to follow fully in 
time the (complicated) trajectory of a Brownian particle. Instead, he introduces the concept 
of a coarse-grained description, defined by a time scale r such that different trajectories 
separated by a time r can be considered independent. No attempt is made to characterize 
the dynamics at a time scale smaller than this coarse-grain time r. A second concept, 
probabilistic in nature, introduced by Einstein is the probability density function (pdf, for 
short), /(A), for the distance A travelled by the Brownian particle in a time interval r. 
/(A) is defined such that f{A)dA is the probability of having a change in position x in the 
interval (A, A + dA). From the fact that /(A) is a pdf it follows the following properties: 

^ (2.1) 
//(A)dA = 1 ^ ' 

One could invoke the law of large numbers to predict a Gaussian form for /(A). However, 
this is not necessary and one only needs to assume the symmetry condition: 

/(-A) = /(A) (2.2) 

We consider an ensemble of Brownian particles. This is characterized by the particle 
number density n(x, t), which is such that n(x, t)dx is the number of particles in the volume 
(x, X + dx) at time t. Prom the assumption that the trajectories separated a time interval r 
are independent, it follows that the number of particles at location x at time t + t will be 
given by the number of particles at location x — A at time t, multiplied by the probability 
that the particle jumps from x — A to x, which is nothing but /(A), and integrated for all 
the possible A values: 

n{x,t + T) = J n{x- A,t)f{A)dA (2.3) 

This is the basic evolution equation for the number density n{x,t). By Taylor expanding 
the above expression, using of the symmetry relation eq.(^]^) and keeping only the lowest 
non-vanishing order terms, one gets the diffusion equation: 

^ = (2.4) 

dt ^ ' 

where the diffusion constant D is given in terms of the second moment of the pdf /(A) by: 

D = ^J A'f{A)dA = ^ (2.5) 

If the initial condition is that all particles are located at x = 0, n{x,t = 0) = N5{x), the 
solution of the diffusion equation is: 

from where it follows that the average position of the Brownian particle is and that the 
average square position increases linearly with time, namely: 



(x(t)) = 
(x(t)2) = 6Dt 



2v _ «n. (2-7) 
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These predictions have been successfuhy confirmed in experiments and contributed to the 
acceptance of the atomic theory. The above results are characteristic of stochastic diffusion 
processes as the ones we wih encounter again in other sections (c.f. Sect. 5). 

Even though Einstein's approach was very successful, one has to admit that it is very 
phenomenological and can not yield, for instance, an explicit expression that allows the 
calculation of the diffusion constant in terms of microscopic quantities. Langevin (1908) 
initiated a different treatment which can be considered in some way complementary of the 
previous one. In his approach, Langevin focused on the trajectory of a single Brownian 
particle and wrote down Newton's equation Force=mass x acceleration. Of course, he 
knew that the trajectory of the Brownian particle is highly erratic and that would demand 
a peculiar kind of force. Langevin considered two types of forces acting on the Brownian 
particle: usual friction forces that, according to Stokes law, would be proportional to the 
velocity, and a sort of "fluctuating" force ^(t) which represents the "erratic" part of the 
force coming from the action of the fluid molecules on the Brownian particle. The equation 
of motion becomes then: 

dv ^ ^ 

m— = —GTrrjav + 4 (2-8) 

T] is the viscosity coefficient and a is the radius of the (assumed spherical ) Brownian particle. 
Langevin made two assumptions about the fluctuating force £,{t): that it has mean and 
that it is uncorrelated to the actual position of the Brownian particle: 

m) = 

{x{t) . m) = im ■ m) = o ^ ■ ^ 

Multiplying eq.(2.8) by x, taking averages with respect to all realizations of the random 
force (,{t), and using the previous conditions on ^(t) one gets: 

Langevin assumed that we are in the regime in which thermal equilibrium between the 
Brownian particle and the surrounding fluid has been reached. In particular, this implies 
that, according to the equipartition theorem, the average kinetic energy of the Brownian 
particle is {mv ^/2) = 2>kBT/2 {ks is Boltzmann's constant and T is the fluid temperature). 
One can now solve very easily eq. (|2.10 ) to find that the particle does not move on the average 



and that, after some transient time, the asymptotic average square position is given by: 

This is nothing but Einstein's diffusion law, but we have now an explicit expression for the 
diffusion coefficient: 

D = (2.12) 

Girrja 

Langevin's random force ^(t) and the Brownian particle position x{t) are examples of 
stochastic processes. It is now time we provide a more precise definition of what a stochastic 



process is. It should be clear that the natural scenario is that of probability theory [15, 16 1. 
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Let us consider a probabilistic experiment E = {S,T,P), where S is a set of possible 
results, is the a-algebra of events, i.e. those subsets of S that are assigned a probability, 
and the real function P : ^ — > is a ir-additive probability. In many occasions, we are 
interested in a real number associated with the experiment result. We call this a random 
variable and denote it by x. In other words, a random variable is an application of the set 
of results into the set of real numbers]^ 

u x[u] 

In many occasions, the outcome u of the experiment is itself a real number, and we simply 
define x(u) = u. In those cases and by abuse of language, the experiment result u is also 
called a random variable. The probability density function f{x) of the random variable x is 
defined such that f{x)dx is the probability that x takes variables in the interval {x,x + dx)^ 
namely the probability of the set {u G S* | x < x[u\ < x -\- dx}. 

We now define a stochastic process x{t) as a family of random variables depending on 
some continuous real parameter t. It is, then, a family of applications: 

x{t) :S^R 

u x[t,u] 

Alternatively, for each experiment result u we might now think of x[t,u] = Xu{t) as a 
function of the parameter t. This is usually the way one considers a stochastic process 
x{t): as a collection of functions Xu{t) each one depending of the outcome u. In most 
applications, t is a physical time and Xu{t) is a trajectory that depends on the outcome 
of some probabilistic experiment S. In this way, the trajectory itself Xu{t) acquires a 
probabilistic nature. 

Arguably, the most well-known example of a stochastic process is that of the random walk 
|l7| . The experiment 5" is a series of binary results representing, for instance, the outcome 
of repeatedly tossing a coin: u = (0, 0, 0, 0, 1, 1, 0, 1, 1, 1,0, 1, 0, 1, 0, 1, . . .) (1 means "heads", 
means "tails"). To this outcome we associate a 1-dimensional real function Xuit) which 
starts at Xu{0) = and that moves to the left (right) at time kr an amount a if the fc-th 
result of the tossed coin was (1). Fig. ( |2.lD shows the "erratic" trajectory for the above 
result u. 

What does one mean by characterizing a stochastic process? Since it is nothing but a 
continuous family of random variables, a stochastic process will be completely characterized 
when we know the joint probability density function for the set {x{ti), x{t2), ■ ■ ■ ,x{tm)}, 
i.e. when we know the function /(xi, . . . , Xm', ti, . . . , tm) for arbitrary m. This function is 
such that 

/(Xi, ...,Xm\ti,.. .,tm)dxi . ..dXm (2.15) 

represents the probability that the random variable x{ti) takes values in the interval 
{xi,xi + dxi), the random variable x{t2) takes values in the interval {x2,X2 + dx2), etc. 
In a different language, we can say that a complete characterization of the trajectory is 
obtained by giving the functional probability density function f{[x{t)]). One has to realize 

^This application must satisfy some additional properties, in particular that the set {u G S | x[u] < x} 
belongs to the a-algebra JF, Va;. 
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(2.13) 



(2.14) 




Figure 2.1: Random walk trajectory Xu{t) associated to the result u of a binary experiment 
as discussed in the main text. 

that a complete characterization of a stochastic process implies the knowledge of a function 
of an arbitrary number of parameters and is very difficult to carry out in practice. In many 
occasions one is happy if one can find simply the one-time pdf f{x;t) and the two-times 
pdf f{xi,X2;ti,t2)- In terms of those, it is possible to compute trajectory averages: 

/oo 
dxx^'fix-t) (2.16) 
-oo 

and time correlations: 

/oo 
dxi dx2 xiX2fixi,X2;ti,t2) (2.17) 
-oo 

It is important to understand that the averages (. . .) are taken with respect to all the 
possible realizations of the stochastic process x{t), i.e. with respect to all the possible 
outcomes u of the experiment. Every outcome u gives rise to a different trajectory x„(i). 
The different trajectories are usually called "realizations" of the stochastic process x{t). 
An alternative way of understanding the previous averages is by performing the experiment 
a (preferably large) number of times M to obtain the results Ui, i = 1,...,M, and the 
different trajectory realizations (t). The averages can then be performed by averaging 
over the different trajectories as: 

, M 

{m = Mj:^uM (2-18) 

1=1 
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and similar expressions for other averages. 

In two very interesting cases does the knowledge of f{x;t) and /(xi, X2; ii, ^2) imply 
the knowledge of the complete pdf f{xi, . . . , x^] ti, . . . , tm) for arbitrary m: (i) Complete 
time independence and (ii) Markov proces^. In a complete time-independent process, the 
random variables at different times are independent and we are able to write: 

f{xi, . . . ,Xm;ti,. . . ,tm) = f{xi;ti)f{x2; ^2) • • • f{Xm; tm) (2.19) 

In the case of a so-called Markov process, the rather general conditional probability 

f ( ■ + \ ■ f f = /(j^l , ■ ■ ■ , Xm] ti, . . . , tm) r)p,\ 

J {Xl , . . . , Xm—1] ) • • • ) t-m—l ) 

is equal to the two-times conditional probability 

ff ,1 , N _ fjXm—l-i Xm'i tm—lj tm) n 

J\Xm7''Tn\Xm—li'^m--l) = 77 |~T \ [■^■■^^J 

J \Xm~ 1 ) '■m— 1 j 

for all times tm > tm-i > ■ ■ ■ > ti. Loosely speaking, the Markov property means that the 
probability of a future event depends only on the present state of the system and not on 
the way it reached its present situation. In this case one can compute the m-times pdf as: 

/(xi, . . . , Xm^ ti, . . . , tm) /2 22") 

f(,Xmj tm\Xm—l: tm—l)f{Xm—l ; tm-l\Xm-2;tm-2) ■ ■ ■ f {X2\ ^2 1^1 ; tl )/(xi ; ti) 

The random walk is an example of a Markov process, since the probability of having a 
particular value of the position at time {k + l)r depends only on the particle location at 
time kr and not on the way it got to this location. 



Another important category is that of Gaussian processes |18| for which there is an 
explicit form for the m-times pdf, namely: 



f{xi, . . . , Xm'i tl, . . . 1 tm) — Y 



m 

-— y ^ (xi — bi)Sij{xj — bj 



(2.23) 



where the parameters of this expression can be related to mean values of the random 
variables as: 

bi = {x{U)) 

{S-% = {x{t,)x{t,)) - {x{t,)){x{t,)) ^^-^^^ 

As a consequence, one very rarely writes out the above form for the pdf, and rather char- 
acterizes the Gaussian process by giving the mean value {x{t)) and the correlation function 
{x{t)x{t')) . Prom the many properties valid for Gaussian processes, we mention that a 
linear combination of Gaussian processes is also a Gaussian process. 

^In fact, complete independence is a particularly simple case of Markov processes. 
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2.2 Stochastic Differential Equations 

A stochastic differential equation (SDE) is a differential equation which contains a stochastic 
process ^{t): 

^ = G(Mt),t,m) (2.25) 

Let us explain a little further what is meant by the previous notation^. G is a given 3- 
variable real function, ^(t) is a stochastic process: a family of functions ^u{t) depending on 
the outcome u of some experiment S. As a consequence a SDE is not a single differential 
equation but rather a family of ordinary differential equations, a different one for each 
outcome u of the experiment S: 

^^ = G{xu{t),t,Ut)) (2.26) 

Therefore, the family of solutions Xu{t) of these differential equations, for different outcomes 
u, constitute a stochastic process x{t). We can say that for each realization ^uit) of the 
stochastic process ^, corresponds a realization Xu{t) of the stochastic process x. The solution 
X becomes then a functional of the process ^. To "solve" a SDE means to characterize 
completely the stochastic process x{t), i.e. to give the m-times pdf /(xi, . . . , Xm', ti, . . . , tm)- 
Again, this is in general a rather difficult task and sometimes one focuses only on the 
evolution of the moments (x(t)") and the correlation function (x(ti)x(t2))- 

When the stochastic process (,{t) appears linearly one talks about a Langevin equation. 
Its general form being: 

dx 

— = q{x,t) + g{x,t)m (2.27) 

(from now on, and to simplify notation, the "hats" will be dropped from the stochastic 
processes). In this case, ^(i) is usually referred to as the "noise" term. A word whose 
origin comes from the random "noise" one can actually hear in electric circuits. Still 
another notation concept: if the function g{x,t) is constant, one talks about additive noise, 
otherwise, the noise is said to be multiplicative. Finally, q{x, t) is usually referred to as the 
"drift" term, whereas g{x, t) is the "diffusion" term (a notation which is more commonly 
used in the context of the Fokker-Planck equation, see later). 

Of course, we have already encountered a Langevin SDE, this is nothing but Langevin's 
equation for the Brownian particle, equation (^). In this example, the stochastic noise was 
the random force acting upon the Brownian particle. The "experiment" that gives rise to 
a different force is the particular position and velocities of the fluid molecules surrounding 
the Brownian particle. The movements of these particles are so erratic and unpredictable 
that we assume a probabilistic description of their effects upon the Brownian particle. 

We will now characterize the process ^(t) that appears in the Langevin equation for the 
Brownian particle. We will be less rigorous here in our approach and, in fact, we will be 
nothing but rigorous at the end! but still we hope that we can give a manageable definition 



*The fact that this equation is a first-order differential equation in no means represents a limitation. If x 
and G are vector functions, this equation can represent any high-order differential equation. For simplicity, 
however, we will consider only first-order stochastic differential equations. Notice that Langevin equation 
for the position of the Brownian particle is indeed a second order stochastic differential equation 
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of the stochastic force ^(t). We first need to define the stochastic Wiener process W{t). 
This is obtained as a suitable hmit of the random walk process [|l6|. The probability that 
the walker is at location x = r a after time t = nr can be expressed in terms of the binomial 
distribution: 

P(x(nr) = ra) = f ^ V"" (2-28) 



From where it follows: 



(x(nr)) = 
(x(nr)^) = na 



2v _ (2.29) 



For n ^> 1 we can use the asymptotic result (de Moivre-Laplace theorem) that states that 
the binomial distribution can be replaced by a Gaussian distribution: 

P{x{nT) <ra) = -+ erf f ^) (2.30) 



2 Vv^. 

(erf(x) is the error function |jl^). We now take the continuum limit defined by: 

n oo, r ^ nr = t 

r — > oo, a — > ra = s (2-31) 
aVr = 1 

with finite t and s. In this limit the random walk process is called the Wiener process W{t) 
and equation ( |2.30| ) tends to: 

P{W{t)<s) = ^ + eri^^^ (2.32) 

which is the probability distribution function of a Gaussian variable of zero mean and 
variance t. The corresponding one-time probability density function for the Wiener process 
is: ^ 

/(T^;i) = -l=exp(-^) (2.33) 



/2^ 2t 

The Wiener process is a Markovian (since the random walk is Markovian) Gaussian process. 
As every Gaussian process it can be fully characterized by giving the one-time mean value 
and the two-times correlation function. These are easily computed as: 

{Wit)) = 

(H^(t)2) = t (2.34) 
{W{h)W{t2)) = min(ti,t2) 



A typical realization of the Wiener process is given in Fig. (p.2|). The Wiener process is 
continuous but it does not have first derivative. In fact it is a fractal of Hausdorff dimension 
2 |20|. 

We will define now the white-noise random process as the derivative of the Wiener 
process. Since we just said that the Wiener process does not have a derivative, it is not 
surprising that the resulting function is a rather peculiar function. The trick is to perform 
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Figure 2.2: A typical realization of the Wiener process done by generating a random walk 
with a large number of steps. 



the derivative before the continuum limit ( |2.3lD is taken. If x{t) is the random walk process, 
we define the stochastic process We{t) as: 

w,{t) = + - ^(^^ (2.35) 
e 

Wf:{t) is a Gaussian process since it is a linear combination of Gaussian processes. Therefore, 
it is sufficiently defined by its mean value and correlations: 

= 

if ti - t2 < -e 

hn(t\.n(f\\ - } aV(^e)(l + (ii-i2)A) if-e<ti-t2<0 (2.36) 

if ii - t2 > e 

which is best understood by the plot in Fig. (^^). If we let now e — > the process 
w{t) = \\m.f:^QWe{t) becomes the derivative of the random walk process. The correlation 
function ( 2.36| ) becomes a delta function (a^/r)5(ti — 12)- If we take now the limit defined 



in eqs.( 2.31 ) the random walk process tends to the Wiener process W{t) and the derivative 
process w{t) tends to S,wit)- the white-noise process. Intuitively, the white-noise represents 
a series of independent pulses acting on a very small time scale but of high intensity, such 
that their effect is finite. The white noise can be considered the ideal limit of a physical 
stochastic process in the limit of very small correlation time r. Formally, the white noise 
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Tea <w,(t)w,(t') > 

1.5 r 




/ n n I \ 

-2-10 1 2 

e"^(t-t') 



Figure 2.3: Correlation function for the derivative of the random walk process, equation 

(ill). 



is defined as a Markovian, Gaussian process of mean value and correlations given by: 

and can be considered the derivative of the Wiener process: 

Ut) = ^ (2.38) 

All this lack of mathematical rigor leads to some problems of interpretation. For in- 
stance, when in a SDE the white-noise appears multiplicatively, the resulting process x{t) 
will be, in general, a non continuous function of time. When this happens, there is an 
ambiguity in some mathematical expressions. Giving a sense to those a priori undefined 
expressions constitutes a matter of sheer definition. The most widely used interpretations 
are those of Ito and Stratonovich |l|, |2|. To make a long story short, we can summarize 
both interpretations as follows: when in a calculation we are faced with an integral 

ds gixis))C^is) (2.39) 



to be computed in the limit h ^ 0, Ito interprets this as: 



gixit)) ds Us) = gixmW{t + h)- W{t)] (2.40) 
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and Stratonovich as: 



+ p = + + ^)) [^(, + _ Wit)] (2.41) 

2 Jt 2 

Although there was much arguing in the past to which is the "correct" interpretation, 
it is clear now that it is just a matter of convention. This is to say, the Langevin SDE 
( |2.27] ) is not completely defined unless we define what do we interpret when we encounter 
expressions such as eq. In some sense, it is related to the problem of defining the 



following expression involving the Dirac-delta function: 

) 

dt 6{t) (2.42) 



This can be deEned as equal to 1 (equivalent to the Ito rule) or to 1/2 (Stratonovich). Both 
integration rules give different answers and one should specify from the beginning which 
is the interpretation one is using. The Stratonovich rule turns out to be more "natural" 
for physical problems because it is the one that commutes with the limit ( 2.3l| ). Unless 



otherwise stated, we will follow in general the Stratonovich interpretation. Moreover, the 
Stratonovich interpretation allows us to use the familiar rules of calculus, such as change 
of variables in an integration, etc. The use of the Ito rules leads to relations in which some 
of the familiar expressions of ordinary calculus are not valid and one needs to be a little bit 
more cautious when not used to it. A consequence of the Ito interpretation is the simple 
result {x{t)S^yj{t)) = which is not valid in the Stratonovich interpretation (see later Eq. 
( |2.56| )). However, there is a simple relation between the results one obtains in the two 
interpretations. The rule is that the SDE 

^ = q(x) + gix)Ut) (2.43) 
in the Ito sense, is equivalent to the SDE: 

^ = qix) - ^gix)g'ix) + gix)Ut) (2-44) 

in the Stratonovich sense. This rule allows to easily translate the results one obtains in 
both interpretations. Both interpretations coincide for additive noise. 

The white noise process is nothing but a physical idealization that leads to some mathe- 
matical simplifications. For instance, it can be proven that the solution x(t) of the Langevin 
equation ( 2.43D is a Markov process if ^(t) is a white-noise process. In any physical pro- 



cess, however, there will be a finite correlation time r for the noise variables. A widely used 
process that incorporates the concept of a finite correlation time is the Ornstein-Uhlenbeck 
noise, ^ou [0- This is formally defined as a Gaussian Markov process characterized by: 

iCouit)) = 

iMmouit')) = ^e-l*-*'!/- ('-^'^ 
Several comments are in order: 

(i) The OU-noise has a non zero correlation time r meaning that values of the noise at 
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different times are not independent random variables. This is a more faithful representation 
of physical reality than the white-noise limit. 

(ii) In the limit r — > the white noise limit (,w{t) (with the correlations given in eq.( |2.37| )) 
is recovered and the corresponding SDE's are to be interpreted in the Stratonovich sense. 

(iii) The OU~noise is, up to a change of variables, the only Gaussian, Markov, stationary 
process. 

(iv) The OU-noise is the solution of the SDE: 

^ = -kou{t) + -Ut) (2.46) 
at T T 

with the initial condition that $,ou{^) is a Gaussian random variable of mean and variance 
(2r)-i. 

2.3 The Fokker-Planck Equation 

We will find now an equation for the one-time pdf for a stochastic process which arises as 
a solution of a SDE with Gaussian white noise. This is called the Fokker-Planck equation 
[^]. It is the equivalent of Einstein's description which focuses on probabilities rather than 
in trajectories as in the Langevin approach. More precisely, we want to find an equation 
for the one-time pdf /(x, t) of a stochastic process governed by the SDE of the Langevin 
form: 

nnf 

- = q{x) + g{x)Ut) (2.47) 

with the white-noise defined in eq. ( [2.371 ). This is to be understood in the Stratonovich 
interpretation. To find an equation for the probability density function f{x;t) we will rely 
upon functional methods |2^. Let us consider first the corresponding deterministic 
initial value problem: 

xit = 0) = xo ^^-^^^ 

The solution of this equation is a (deterministic) function x{t,xo). We can think of x{t) as 
a random variable whose probability density function p(x; t) is a delta-function: 

p{x;t) = 6{x - x{t,xo)) (2.49) 

We can now turn x{t) very easily into a stochastic process by simply letting the initial 
condition xq become a random variable. For each possible value of xq we have a different 
solution x(t, Xo), i.e. a random process. In this case the probability density function p{x,t) 
is obtained by averaging the above pdf over the distribution of initial conditions: 

p{x;t) = {6{x-x{t,xo)))^^ (2.50) 

But it is well known from mechanics that the density p{x; t) satisfies Liouville's continuity 
equation [p^ : 

| + |:(xp)=0 (2.51) 
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or, using, eq.(p.48|): 



dp d 



(2.52) 



If we consider now the full SDE ( ^.47 ) we can repeat the above argument for a given 
realization of the noise term. The probability density function f{x;t) will be the average 
of p{x; t) with respect to the noise distribution: 



fix; t) = {p{x; = {d{x - x{t, xo)))xo,€» 
where p satisfies the Liouville equation ( |2.51| ) and, after substitution of ( p.47| ): 

By taking averages over the noise term we get: 



dl 
dt 

The averages are done by using Novikov's theorem [p4| : 

{mpmm) 







(2.53) 



(2.54) 



(2.55) 



(2.56) 



which applies to any functional p{[^{t)\) of a Gaussian process ^(t) of zero mean, (^(t)) = 0. 
For the white-noise process: 



{ut)p{[um) 



1 



5p 



2 \5iu,{s) 

By using the functional calculus, this can be computed as: 
Sp \ I 5x{t) 



s=t 



s=t 



s=t 



Sp \ _ d / 6x{t) 



P) 



s=t 



By using the formal solution of eq. ( 2.47| ) : 



we get: 



and 



x{t) = XQ + 1^ dsq{x{s))+ I ds g{x{s))^w{s) 
6x(t) 



9{x{t)) 



s=t 



dp 



-^b(x)(.)] = -^[.(.)/] 



(2.57) 
(2.58) 

(2.59) 

(2.60) 
(2.61) 



By substitution in ( 2.57| ) and ( |2.55| ) we get, finally, the Fokker-Planck equation for the 
probability density function: 



dl 
dt 



(2.62) 
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The above functional method to derive the Fokker-Planck equation from the Langevin 
equation is very powerful and can be extended to other situations such as the multivariate 
case and SDE with colored (Ornstein-Uhlenbeck) noise. For colored noise the problem 
is nonmarkovian and no exact Fokker-Planck equation can be derived for the probability 
density. Still one can obtain different approximate equations for the probability density 
and for the two-time probability density [p5|]. 

Much work has been devoted to finding solutions to the Fokker-Planck equation. For 
the stationary equation, dP/dt = 0, the solution can be always reduced to quadratures in 
the one-variable case. The multivariate equation is more complicated and some aspects are 
discussed in Sect. 4.2. 



2.4 Numerical generation of trajectories 

In the Langevin's approach to stochastic processes, a great deal of relevance is given to the 
trajectories. It is of great importance, therefore, when given a SDE of the form of eq.(|2^). 



to be able to generate representative trajectories. This is to say: to generate the functions 
Xu{t) for different outcomes u of the experiment. If we generate, say M results Uj, the 
averages could be obtained by performing explicitly the ensemble average as indicated in 



equation (2.1 



We now explain the basic algorithms to generate trajectories starting from a SDE of the 



Langevin form [21, 26, 27, 



2.4.1 The white noise case: Basic algorithms 

Euler's algorithm is by far the simplest one can devise to generate trajectories. We will 
explain it by using a simple example. Let us consider the SDE: 

m = f{t)+iu,{t) (2.63) 

If fact, this equation is so simple that it can be solved exactly. For a given realization of 
the noise term, the solution is: 

x{t) = x(0) + /* f{s) ds + /* C^(s) ds = x{Q) + F{t) + W{t) (2.64) 

JO JO 

W{t) is the Wiener process. Hence, we conclude that the stochastic process x{t) is Gaussian 
and is completely characterized by: 

(x(t)) = x(0) + F(t) 
{x{t)x{t')) = (x(0) +F(t))(x(0) +F(t')) +min(t,t') ^^^^> 

However, let us forget for the moment being that we can solve the SDE and focus on a 
numerical solution in which we generate trajectories. We do so by obtaining x{t) at discrete 
time intervals: 

x{t + h) = x{t) + ^l"^^ x{s) ds 

= xit) + J,'+''f{s)ds + Jl+''Us)ds (2.66) 

= x{t) + fh{t) + Wh{t) 
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here, Wh{t) = W{t + h) — W{t) is the difference of the Wiener process at two different times 
and it is, therefore, a Gaussian process. Wh{t) can be characterized by giving the mean and 
correlations: ^ ^ 

{whit)) = [ {^^{s))ds = (2.67) 



{wh{t)whit')) = Jl""^ ds jf,""^ du{iu,{s)iw{u)) . ^ 

= h dsj^, dud{s-u) 

The integral is an easy exercise on 6 function integration: we can assume, without loss 
of generality that t' > t. li t' > t + h the integral is since there is no overlap in the 
integration intervals and the delta function vanishes. If t < t' < t + h the double integral 
equals the length of the overlap interval: 

ft+h 

{whit)wh{t')) = ds = t-t' + h (2.69) 



In particular, we notice the relation 

{whitf) = h (2.70) 

It is important to realize that, if tj = i /i, tj = j h are the times appearing in the recurrence 
relation eq.( |2.66| ), we have: 

{wh{ti)wh{tj)) = h6ij (2.71) 

We introduce now a set of independent Gaussian random variables u{t) defined only for the 
discrete set of recurrence times, t = 0, h, 2h, . . ., of mean zero and variance one: 

{uit))=0 , {u{t^)) = l 
{u{t)u{t')) = , t^t' ^ ' ' 

There is a vast amount of literature devoted to the question of generation of random num- 
bers with a given distribution |2£, 31, 32|. The set of numbers u{t) can be generated by 
any of the standard methods available. One of the most widely usedQ is the Box-Muller- 
Wiener algorithm: if ri and r2 are random numbers uniformly distributed in the interval 
(0, 1) the transformation: 

gi = v /-21n(ri) cos(27rr2) . . 

52 = \/-21n(ri)sin(27rr2) ^ ' ' 

returns for gi, g2 two Gaussian distributed random numbers of mean zero and variance 
one. This, or other appropriate algorithm, can be used to generated the set of Gaussian 
variables u{t). In terms of this set of variables we can write: 

Wh{t) = h^^^ u{t) (2.74) 

Finally, the recurrence relation that generates trajectories of eq.( ^.63| ) is the Euler algo- 
rithm: 

x{t = 0) = xo 

x{t + h) = x{t) + hit) + h^^ u{t) ^^-'^^ 

^Although not the most efficient. See ref. jsij for a comparison of the timing of the different algorithms 
and the description of a particularly efficient one. 
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For the deterministic contribution we can approximate fhit) ~ hf{t) from where it fohows 
that the deterministic contribution is of order and successive contributions go as 
etc. On the other hand, the contribution coming from the white noise term is of order /i^/^ 
and, in general, successive contribution will scale as h^, h^^"^, etc. 

With the experience we have got by solving the previous simple example, let us now 
tackle a more complicated case. Let us consider the following SDE: 

x{t) = q{x) + g{x)^y,{t) (2.76) 

At this step, one might ask why, given that for a particular realization of the noise term 
the equation becomes an ordinary differential equation (ode), we need special methods to 
deal with this kind of equations. The answer lies in the fact that all the methods developed 
to deal with ODE's assume that the functions appearing in the equation have some degree 
of "well-behaveness" . For instance, they are differentiable to some required order. This 
is simply not the case for the white-noise functions. Even for a single realization of the 
white-noise term, the function is highly irregular, not differentiable and, in our non rigorous 
treatment, is nothing but a series of delta-functions spread all over the real axis. This is 
the only reason why we can not use the well known predictor-corrector, Runge-Kutta and 
all the like methods without suitable modifications. If our SDE happened to have smooth 
functions as random processes we could certainly implement all these wonderful methods 
and use all the standard and very effective routines available! However, this is usually not 
the case: the stochastic contribution is a non analytical function and we must resource to 
new and generally more complicated algorithms. The answer lies in integral algorithms: 
whereas the derivatives of the white-noise function are not well defined, the integrals are 
(the first integral is the Wiener process, which is a continuous function). We look for a 
recursion relation by integration of the SDE ( 2.7(j ): 



rt+h j-t+h 

x{t + h)-x{t)= / q{x{s))ds+ / g{x{s))iui{s) ds (2.77) 



t Jt 



Now we assume that the functions q{x) and g{x) are differentiable functions and we Taylor- 
expand q{x{s)) and g{x{s)) around x = x(t): 



q{x{s)) = q{x{t)) + ^ 

g{x{s)) = g{x{t)) + ^ 
ax 



{x{s) - x{t)) + 0[[x{s) - x{t)y] (2.78) 

x(t) 

{x{s) - x{t)) + 0[{x{s) - x{t)f] (2.79) 



x{t) 

Substitution of the lowest possible order of these expansions, i.e., q{x{s)) = q{x{t)), g{x{s)) 



g{x{t)) in eq.( 2.77 ) yields: 

x{t + h)- x(t) = hq{x{t)) + hO[x{s) - x{t)] + Wh{t)g{x{t)) + Wh{t)0[x{s) - x{t)] (2.80) 
where, as in the simple previous example, 

Wh{t) = ds C^(s) = h^/^ u{t) (2.81) 
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is of order /i^/^. To the lowest order {h^^'^) we have: 

x{s) - x{t) = g{x{t)) r dv CUv) = 0[h'/^] 



(2.82) 



Therefore we need to go one step further in the Taylor expansion of the function g{x{s)) 
and the next -order contribution to ( p.77|) is: 



g'{x{t)) //+'^ ds {x{s) - x{t))Us) + 0[x{s) - x{t)]^WH{t) 
g'{x{t))g{x{t)) Jl+^ ds // dv Us)Uv) + 0[h^/^] 

The double integral can be easily done by changing the integration order: 



(2.83) 



t+h 



ds / dv £,w{s)C^^{v) 



t+h 



t+h 



dv / ds £,wis)(,w{v) 



t+h 



ds 



t+h 



(2.84) 

where to go from the first expression to the second we have exchanged the integration order 
and, to go from the second to the third, we have exchanged the variables names v ^ s. 
Since the first and third integrals are equal they also equal one half of their sum, and the 
previous integrals can be replaced by: 



I rt+h 



dv 



t+h 



I rt+h 



ds Cwis) 



t+h I 

dv CM = o K(*)] (2.85) 



Putting all the bits together, we arrive at the desired result which is 0: 



x{t + h)= x(t) + h}-/'^g{x{t))u{t) + h 



Qix{t)) + -gixit))g'ixit))uit)' 



This recurrence relation is known in the literature as Milshtein method |21, 33]. If the noise 
is additive: g'{x) = 0, then the resulting algorithm is called the Euler algorithm: 



x{t + h)= x{t) + hq{x{t)) + g{x{t))h^l'^u{t) + 0[h^l'^] 



(2.87) 



Sometimes, though, the name "Euler algorithm" is also given to a modification of the 
Milshtein algorithm in which the u{t)^ term is replaced by its mean value: {u{t)'^) = 1: 



x{t + h)= x{t) + g{x{t))h^l^u{t) + h[q{x{t)) + \g{x{t))g' {x{t))] + 0[/i3/2] 



(2i 



This "Euler algorithm" is the one appearing naturally when one does the numerical inte- 
gration of the SDE in the Ito formalism, although at this level it has to be considered just 
as an approximation, unnecessary, to the Milshtein method. 

In the previous expressions, the correction to the recurrence algorithm is said to be of 
order 0[/i^/'^]. Let us explain a little further what is meant by this expression. We use the 
following notation: we call x{t) the values obtained from the numerical integration following 
the Milshtein method: 



x{t + h)= x{t) + h^/'^g{x{t))u{t) + h 



1 



(2.89) 



®It is instructive to realize that the same result can be obtained by use of the Stratonovich rule Eq.( |2.4 
in equation (2.77) 
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and we want to compare x{t + h) with the exact value x{t + h) which is obtained by exact 
integration of the differential equation starting from x{t). What we have proven is that the 
mean-square error (averaged over noise realizations) for the trajectories starting at x{t) is 
of order 0[h^]: 

{{x{t + h) - xit + h)f) = 0[h^] (2.90) 

One says that the Milshtein algorithm has a convergence for the trajectories in mean square 
of order h'^. This is related, but not the same, as the order of convergence of the n-th order 
moment, Z?„(/i), which is defined as: 

D„(/i) = (x(t + /i)") - (x(t + /i)") (2.91) 

the averages are done starting from a given x{t) and averaging over noise realizations. For 
the Milshtein algorithm one can prove p7| : 

Dn{h) = 0[h'^] (2.92) 

Which means that, when computing moments, the Milshtein algorithm makes an error of 
order h'^ in every integration step. Of course, in a finite integration from t = to a time 
t = kh, the total error will be multiplied by the number of integration steps, k, which gives a 
contribution fcO[/i^] = {t/h)0[h'^] = 0[h]. And we can write the following relation between 
the exact value a;(f)" and the value x(t)" obtained when using the Milshtein approximation 
starting from an initial value x(0): 

{x{t)'') = {x{tr) + 0[h] (2.93) 

In practice, what one does (or rather, what one should do!) is to repeat the numerical 
integration for several time steps hi, /12, ... and extrapolate the results towards /i = by 
using a linear relation {x{t)"') = (x(t)") + a • /i 

The question now is whether we can develop more precise algorithms while preserving 
the structure of the Milshtein method, i.e. something like: 

x{t + h) = x{t) + CiWhit) + C2h + C:iWh{tf + Cihwh{t) + C^Whitf + ... (2.94) 

The (negative) answer was given by Riimelin |^^, who stated that higher order algorithms 
necessarily will imply more random processes, say etc. But the problem lies on 

the fact that those random processes are not Gaussian and have non-zero correlations being, 
in general, very difficult to generate accurately. As a conclusion, the Milshtein algorithm 
appears as the simplest alternative for integrating a SDE. However, and as a recipe for 
the practitioner, Runge-Kutta type methods offer some advantages at a small cost in the 
programming side. These methods will be explained in a later section. As a final remark, 
the Milshtein method can also be used for the SDE: 

x{t)=q{t,x)+g{t,x)^^{t) (2.95) 

in which the diffusion an drift terms depend explicitly on time t. 
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2.4.2 The Ornstein Uhlenbeck noise 

We turn now to the numerical generation of trajectories for a SDE with colored noise, in 
particular of the Ornstein-Uhlenbeck form as defined in equation (2.45). First, we explain 
how to generate realizations of the OU process itself, ^ouit)^ and later we will see their use 
in SDE's. 

Equation (2.46) can be solved exactly (it is a linear equation) and the solution actually 
tells us how to generate trajectories of the OU-process. The solution is: 



Couit + h) = ^ou{t)e-''/^ + Hh{t) 
where we have introduced the random process Hh(t) as: 



Hh{t) = T-'e 



t+h 



ds iw{s)e: 



s/t 



(2.96) 



(2.97) 



Using the white noise properties, eq. (|2.37| ), it is an easy exercise to prove that Hh{t) is a 
Gaussian process of mean value zero: 



and correlations: 

{Hh{t)Hh{t')) 



(2r)-i 




{Hh{t)) = 



-\t-t'\/T _J\t~t'\-2h)/r 



if \t -t'\<h 
if \t -t'\> h 



(2.98) 



(2.99) 



The important part is to realize that, for the times ti = ih that appear in the recurrence 
relation eq.( |2.96| ), the correlations are: 



{Hhik)Hhitj)) = (2r)-i 



1 



-2h/T 



(2.100) 



If we introduced a set of independent Gaussian variables u{t) of zero mean and variance 
unity, the process Hh{t) can be generated as: 



Ht,{t) = v/(l - e~^h/ry^2T)u{t) 



(2.101) 



And the final recurrence relation to generate trajectories of the Ornstein-Uhlenbeck noise 
is: 

/ ^ouiO) 



Cou(t)e-''/^ + V^(l - e-2Vr)/(2r)n(t + h) 
Let us consider now an SDE with OU noise. Let us start again by the simplest example: 

dx{t) 



dt 



f{t) + ^ou{t) 



(2.103) 



We use an integral algorithm: 

i-t+h 



rt+h rt+h 

x{t + h)= x{t) + / ds f{s) + / ds Couis) = x{t) + fh{t) + Qhit) (2.104) 
Jt Jt 
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The stochastic contribution gh{t) is a Gaussian process characterised by the following mean 
value and correlations: 



{9h{t)) 
{9h{t)gh{t')) 



cosh ( ^ 



1 



exp 



(2.105) 



(valid for all the times t, t' appearing in the recursion relation eq.( p.l04 ).) To order 0[/i^], 
the correlations become: 



{ghit)ghit')) 



2t 



exp 



\t-t'\ 



(2.106) 



from where it follows that the process ghit) is nothing but h times an Ornstein-Uhlenbeck 
process: gh{t) = h^ou{t)- Summarizing, to order /i^ the algorithm to integrate numerically 
the SDE eq. (|2l03|) is: 



x{t + h)= x{t) + fh{t) + hiou{t) + 0[h? 



(2.107) 



where S,ou{t) is generated by the use of eqs. (12.1021) . 

If one needs more precision in the stochastic part, one step further in the integration of 
the equation can be achieved by generating exactly the process gh{i) We define the 

process: 



G{t) 



in terms of which: 



gh{t) = G{t + h) - G{t) 



(2.108) 



(2.109) 



Since dG{t)/dt = £,ou{t) and iou{t) satisfies the differential equation eq. (|2.46| ) we can write 
down the following equation for G{t): 



(fG{t) , 1 dG{t) 
~~dP ^ T dt 

G(0) 

dG(t) I 





eoc/(0) = Co 



(2.110) 



whose solution is: 

G{t) = rCo - T^oe-*/" + /* ds i^{s) - e"*/^ /* ds e^/"C^(s) 

Jo Jo 

From where it follows the recursion relation: 

gh{t + h)= pgh{t) - ph{t) + h{t + h)- f2{t) + f2{t + h) 



(2.111) 



(2.112) 



The initial condition is that gh{0) is a Gaussian variable of mean and variance given by 
equation ( ^.105| ) for t = t' = 0. This can be written as: 



gh{0) 



cosh I — 



u 



(2.113) 
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where u is a Gaussian random number of zero mean and unit variance. In equation ( p.ll2| ) 
we have introduced the following definitions: 



P 

hit) 

f2{t) 



(2.114) 



The processes fi{t) and /2(t) are correlated Gaussian processes, whose properties, for the 
times ti = ih appearing in the recursion relations, are given by: 



ifiiti) 

{f2{ti) 

{f2{ti)f2{tj) 
mi)f2{h) 








(2.115) 



2 '^^J 

r(l - p)6ij 



It is possible to generate the processes fi{t) and f2{t) satisfying the above correlations by 
writing them in terms of two sets of independent Gaussian distributed random numbers of 
zero mean and unit variance, ui{t), U2{t): 



fi{t) = aiui{t) 

h{t) = PlUi{t) + (32U2{t) 



(2.116) 



where the constants ai, (3i and /?2 are chosen in order to satisfy the correlations ( |2.115 ): 

(2.117) 



ai 
Pi 



t(1-p) 
Vh 



P2 



r{l-p) 
2 



2t 



+ P 1 + 



2t 



In summary, the process gh{t) is generated by the recursion relation ( |2.112| ) with the initial 
condition (|2.113 ) and the processes fi{t) and /2(i) obtained from the relations ( ^.116|) . 
We now consider a more general equation: 



dx 

— = q{t,x) + g{t,x)iou{t) 
at 



We start again by an integral recursion relation: 



/•t+Zi rt+n 

x{t + h) = x{t)+ f dsq{s,x{s))+ / ds g{s,x{s))iou{s) 
Jt Jt 

By Taylor expanding functions q and g one can verify that at lowest order: 
x{t + h) = x{t) + hq{t, x{t)) + gh{t)g{t, x{t)) + 0[h?\ 



rt+h 



(2.118) 



(2.119) 



(2.120) 



where gh{t) is the process introduced before. As explained, ghit) can be generated exactly. 
Alternatively, and in order to make the resulting algorithm somewhat simpler, one can 
replace gh{t) by h^ou{t) + 0[h'^] without altering the order of convergence of the algorithm. 
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However, these algorithms suffer from the fact that they do not reproduce adequately the 
white-noise Milshtein method. This is to say: always the integration step h has to be kept 
smaller than the correlation time r. If one is interested in the limit of small values for the 
correlation time r (in particular, if one wants to consider the white noise limit r ^ 0), it is 
better to turn to the Runge-Kutta methods which smoothly extrapolate to the Milshtein 
method without requiring an extremely small time step. 



2.4.3 Runge Kutta type methods 

We focus again in the SDE with white noise: 

x{t)=q{t,x)+g{t,x)i^{t) (2.121) 

We will develop now a method similar to the second-order Runge-Kutta (RK) method 
for solving numerically ordinary differential equations. As we said before, the particular 
features of the white-noise process prevent us from simply taking the standard R-K methods 
and we need to develop new ones. 

Let us recall briefly how a Runge-Kutta method works for an ordinary differential equa- 
tion: 

^ = g(t,x) (2.122) 

Euler method: 

x{t + h)= x{t) + hq{t, x{t)) + 0[h'^] (2.123) 

can be modified as: 

h 

x{t + h) = x{t) + -[q{t,x{t)) + q{t + h,x{t + h))] (2.124) 

Of course, this is now an implicit equation for x{t + h) which appears on both sides of the 
equation. RK methods replace x[t + h) on the right hand side by the predictor given by 
the Euler method, eq.( p.l2^ ) to obtain an algorithm of order 0[/i^]: 

x{t + h)= x{t) + ^[q{t, x{t)) + q{t + h, x{t) + hq{t, x{t)))] + 0[h^] (2.125) 



This is usually written as: 



k = hq{t,x{t)) 
x{t + h) = x{t) + ^[q{t,x{t))+q{t + h,x{t)+k)] 



(2.126) 



The same idea can be applied to the SDE (|2.121|) . Let us modify Euler 's method, eq. (|2^ ) 
(which we know is a bad approximation in the case of multiplicative noise anyway) to: 

xit + h) = x{t) + !l[q{t,x{t))+q{t + h,x{t + h)] + 
!^^[g{t,x{t))+git + h,x{t + h))] 



24 



And now replace x{t + h) on the right-hand-side by the predictor of the Euler method, 
eq.( |2.87| ) again. The resulting algorithm: 



k = hq{t,x(t)) 
I = h^/^u{t)g{t,x{t)) 
x{t + h) = x{t) + ^[q{t,x{t))+q{t + h,x{t)+l + k)]+ ^^'^^^^ 

x{t)) + g{t + h, x{t) +k + l)] 



is known as the Heun method |2^. To study the order of convergence of this method one 
can Taylor expand functions q{t, x) and g{t, x) to see that one reproduces the stochastic 
Milshtein algorithm up to order h. Therefore, from the stochastic point of view, the Heun 
method is of no advantage with respect to the Milshtein method. The real advantage of 
the Heun method is that it treats better the deterministic part (the convergence of the 
deterministic part is of order h^) and, as a consequence, avoids some instabilities typical of 
the Euler method. 



Similar ideas can be applied to the integration of the SDE (2.118) with colored noise 
(111). It can be easily shown that the RK type algorithm: 



k = hq{t,x{t)) 
I = gh{t)g{t,x{t)) 
x{t + h) = x{t) + ^[q{t,x{t))+q{t + h,x{t)+k + l)]+ ^^■^^^> 
SlM[g^t,xit))+g{t + h,x{t)+k + l)] 

correctly reproduces the algorithm ( |2.120[) . Moreover, when the stochastic process gh{t) is 
generated exactly as explained before, this algorithm tends smoothly when r — > to the 
Milshtein algorithm for white noise without requiring an arbitrarily small integration step 
h. 



2.4.4 Numerical solution of Partial Stochastic Differential Equations 

We end this section with the basic algorithms for the generation of trajectories for partial 
stochastic differential equations (PSDE). We will encounter several examples in Sects. 5 
and 6. In general, one has a field A{f,t), function of time t and space r, that satisfies a 
PSDE of the form: 

8A 

— = G[A, VA, V^A, ^(f, t)] (2.130) 

Where G is a given function of the field A and its space derivatives. For the stochastic field 
i^(r, t) usually a white noise approximation is used, i.e. a Gaussian process of zero mean 
and delta-correlated both in time and space: 

mmr',t'))=5{f-r')6{t-t') (2.131) 

The numerical solution of ( |2.130| ) usually proceeds as follows: one discretizes space r ^ ri 
by an appropriate grid of size Ar. The index i runs over the N lattice sites. Usually, but 
not always, one considers a d-dimensional regular lattice of side L, such that N = (L/Ar)*^. 
The elementary volume of a single cell of this lattice is (Ar)'^. Next, one replaces the fields 
A{r,t), (,{f,t) by a discrete set of variables. For the field A{fi,t) we simply replace it by 
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A{fi,t) Ai{t). For the white noise, we have to consider the delta-correlation in space 
and use the substitution: 

m,t) ^ iAr)-^/^Ut) (2.132) 

which comes from the relation between the Dirac 6{f) and the Kronecker 6i functions. 
In this expression ^j(i) are a discrete set of independent stochastic white processes, i.e. 
Gaussian variables of zero mean and correlations: 

mMt')) =Sij6{t-t') (2.133) 

In the simplest algorithms, the field derivatives are replaced by finite differences]]. For 
instance: if the d-dimensional regular lattice is used for the set fi, the Laplacian V'^A[ri,t) 
can be approximated by the lattice Laplacian: 

V^A{f,, t) « (Ar)-2 [Aj{t) - A,{t)\ (2.134) 

Where the sum runs over the set of 2d nearest neighbors of site rj. With these substitutions 
the PSDE ( |2.130D becomes a set of coupled ordinary differential equations: 

= G,(Ai(t), . . . , AM{t);Ut). CN{t)), i = l,...,N (2.135) 

In most occasions, these equations are of the generalized Langevin form: 

dAi{t) 



dt 



(l^m+Y.9^m)m (2-136) 



[A\ denotes the set [A\ = (Ai, . . . , Aat) and gj([^]), gij([A]) are given functions, maybe 
depending explicitly on time. The numerical integration of (|2.136 ) proceeds, as in the 
single variable case, by developing integral algorithms 28|. It turns out, however, that 
in the most general case it is very difficult to accurately generate the necessary stochastic 
variables appearing in the algorithms. That is the reason why one rarely goes beyond 
Euler's modification of Milshtein's method (eq. ( p. 88 ), which reads [p8| : 

Ai{t + h) = Ai{t) + h^'^Y.9^M{t)]>jit) + ^ 



qi{[A{t)\) + - }_^gjk{[A{t)\) 



(2.137) 

Ui{t) are a set of independent random variables defined for the time 0, /i, 2h of the recurrence 
relation, with zero mean and variance one: 

(ui(t))=0 , {ui{t)uj{t)) = 5ij ^9 1^8^ 

(n,(t)n,(t')) =0 , t^t' ^^•'''^^ 

which can be generated by the Box-Muller~Wiener or an alternative algorithm. We stress 
the fact that the functions gjj([A]) are of order (Ar)^'^/^ due to the substitution ( p.l32| ). 

^Alternatively, one could use Fourier methods to compute the Laplacian or other derivatives of the field; 
see, for instance, |37[. 
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For small Ar, this usually demands a small time-integration step for the convergence of 
the solution ||, |4|. 

An important case in which one can use straightforward generalizations of the Milshtein 
and Heun methods is that of diagonal noise, i.e. one in which the noise term does not 
couple different field variables, namely: 



dAjjt) 
dt 



qi{[A]) + g,{Ai)ii{t) 



In this case, the Milshtein method reads: 
Ai{t + h)= Ai{t) + gi{Ai{t))h}'^Ui{t) + h 



u{[Am + -g,{A,{t))9\{Mt)>i{tf 



(2.139) 



(2.140) 



The Heun method is also easily applied in the case of diagonal noise: 

h = hqi{[A{t)]) 
k = h^l-'ui{t)gi{[A{t)]) 
A,{t + h) = A,{t) + ^[qi{[A{t)])+q,{[A{t)+l + k])] + 

^^^%^[5i(A,(t)) + g,{A,{t) +h + k)] 



(2.141) 



2.5 A trivial (?) example: The linear equation with multiplicative noise 

In previous sections, we have shown how it is possible to generate trajectories from a 
given stochastic differential equation and how these trajectories can help us to perform 
the necessary averages. Of course, if one has computed the probability density function 
(by solving, may be numerically, the Fokker-Planck equation |^]) one can also perform the 
same averages. However, there are cases in which much information can be learnt about 
the solution of an stochastic differential equation by looking at individual trajectories. 
Information which is not so easy to extract from the probability density function. To 
illustrate this point we will analyse in some detail the apparently simple SDE: 

— = [a + aUt)]x{t) (2.142) 

with ^w{t) a white noise process with correlation given by eqs.( p.37|) . In fact, this linear 
stochastic differential equation is so simple that it can be solved explicitly: 

x(t) = x(0)e"*+"^(*) (2.143) 

where W{t) is the Wiener process. From the explicit solution and using the known prop- 
erties of the Wiener process we can compute, for instance, the evolution of the mean value 
of x(t): 

{x{t)) = (x(0)e'^*+"^W) = x(0)e("+"'/2)t (2.144) 
where we have used the result 

(e^) = e<^'>/2 (2.145) 

valid for a Gaussian variable z of zero mean. From equation ( p.l44| ) it follows that the 
mean value of x{t) grows towards infinity for a > —a'^/2 and decays to zero for a < —a'^/2. 
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Figure 2.4: Evolution of the mean value for the linear equation ( p.l42| ) for a = —0.1, 
a = 1, xo = 0.01. The solid line is the exact result eq. ( |2.144 ) and the symbols and error 



bars come from a numerical solution of the equation. Notice that although the numerical 
solution follow closely the analytical result, the error bars greatly increase with time. This 
is a clear signature of very strong fluctuations, see the main text. 

Fig. (^^) shows this exponential behaviour of the mean value together with some results 



obtained by numerically integrating the SDE (2.142). This exponential growth is at variance 



to what happens in the deterministic case, o" = 0, for which x{t) decays to for a < and 
grows for a > 0. One would say, in this case, that there has been a shift in the critical 
value of the parameter a to induce a transition from the state x = to the state x ^ 0, 
representing perhaps the transition from a disordered to an ordered state. Of course, it is 
of no concern the fact that {x{t)) tends to infinity. In a real system, they will always be 
saturating terms that will stop growth of x{t) and will saturate {x{t)) to a finite value. For 
instance, a realistic equation could be one with a saturating cubic term: 

^ = [a + aUt)Mt)-<tf (2.146) 

The conclusions of our simple linear analysis would then say that the stationary state of 
the non-linear equation (2.146) is for a < Oc and non-zero for a > a^, being = 



a "critical" value for the control parameter o. Obvious conclusion ... or is it? Well, let us 
have a closer look. 

A first signature that something is missing in the previous analysis is that we could 
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Figure 2.5: Evolution of the mean value for the non-linear equation (2.146) with the same 
values of the parameters than in Fig. (2^). The results, with the error bars, come from a 
numerical simulation of eq.( 2.146| ). 



repeat it as well for the evolution of the n-th moment with the result: 

{x{t)"') = ((x(0)e''*+^W)") = 2;(0)"e"(''+"'^'/2)t (^2.147) 

We see that the n-th moment of the linear equation ( p.l42| ) diverges at a critical value 
0'c,n = —na'^/2 that depends on the order n of the moment. By repeating the arguments 
sketched above, one would conclude that the asymptotic value of the n-th moment of the 
non-linear equation ( |2.146| ) would change from zero to a non-zero value at Oc,n and hence, 
the location of the putative transition from order to disorder would depend on n, which 
does not make much sense. 

The solution of the associated Fokker-Planck equation associated to the non-linear 
Langevin equation ( |2.146| ) and the complete calculation of the moments is given 

in ||4^, ^ w ith the conclusion that all moments tend to for a < 0. This is illustrated 
in Fig. (|2.5| ) in which we plot the time evolution for the mean value {x{t)) for a situation 
in which the linear equation explodes. Since the mathematics used in this proof do not 
provide much physical intuition, it is instructive to rederive the same result by studying 
the behavior of individual trajectories. This will also shed light on the relation between the 
behavior of the moments in the linear and the non-linear equations. 

If we look at the solution, eq.( 2.143| ), of the linear equation, we can say that the determin- 



istic contribution at will always dominate for large time t over the stochastic contribution, 
W{t), which, according to eq.( 2.34| ) is of order t^/^. Hence, for large t, and for a < every 
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trajectory will go to zero, and consequently {x{t)^) = 0, in contradiction with previous re- 
sults, in particular equation ( p.l47| ). If fact, the statement that at dominates over W{t) and 
hence x{t) decays to zero is not very precise since at is a number and W{t) is a stochastic 
process. A precise statement is that the probability that x{t) decays to zero (i.e. that it 
takes values less than any number e) tends to one as time tends to infinity: 

Ve, lim Probability < e) = 1 (2.148) 

Remember that W{t) is a stochastic Gaussian variable and, hence, it has a finite probability 
that it overcomes the deterministic contribution at at any time. What the above relation 
tells us is that, as time tends to infinity, this probability tends to zero, and we can say, 
properly speaking, that every trajectory will tend to zero with probability one. However, for 
any finite time there is always a finite (however small) probability that the stochastic term 
W{t) overcomes the deterministic contribution at by any large amount. This is illustrated 
in Fig. ( |2.6| ) in which we plot some typical trajectories for the linear equation. In these 
figures it is important to look at the vertical scale to notice that, although every trajectory 
decays to zero as time increases, there are very large indeed fluctuations for any given 
trajectory. 

In brief, for the linear equation ( |2.142 ) there are unlikely trajectories (with decreasing 



probability as time increases) that become arbitrarily large. It is the contribution of those 
small probability, large amplitude trajectories, which make the moments diverge. 

Now we can look at the non-linear problem with other eyes. The non-linear term will 



suppress those large fluctuations in the trajectories one by one, not just on the average [44|. 



We conclude, then, that the mean value of any moment {x{t)") will tend to zero for a < 0, 



in agreement, of course, with the exact result in |42, 43 1. We can say that the presence of 



the multiplicative noise makes no shift in the transition from ordered to disordered states. 



Fig. (2/7) shows some trajectories in the non-linear case to see that, effectively, there are 
no large fluctuations as one can see comparing the vertical scale of Figs. (2.6) and (2.7). In 
Sect. 6 we will analyze other situations in which noise can actually induce an order/disorder 



transition. 



3 Transient stochastic dynamics 

In many physical situations noise causes only small fluctuations around a reference state. 
Of course, noise effects are more interesting when this is not the case, as for example when 
there is a mechanism of noise amplification. A typical situation of this type is the decay 



of an unstable state |45]: a given system is forced to switch from one state to another by 
changing a control parameter. After the change of the parameter, the system, previously 
in a stable state, finds itself in a state which is an unstable fixed point of the deterministic 
dynamics, and it is driven away from it by fluctuations. Noise starts the decay process 
and fluctuations are amplified during the transient dynamics. Examples of such situations 
include spinodal decomposition in the dynamics of phase transitions and the process of 
laser switch-on |4^. The latter essentially consists in the amplification of spontaneous 
emission noise. We will discuss laser switch-on in the second part of this section after some 
general methodology is presented. A related interesting situation, which we will not discuss 
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Figure 2.6: Series of five trajectories x{t) for the linear multiplicative problem to show 
that there are large amplitude fluctuations. Same parameters than in Fig. 
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Figure 2.7: Series of five trajectories x{t) for the non-linear multiplicative problem 
eq.( |2.14^ ) to show that the large amplitude fluctuations are suppressed by the non-linear 
terms for each trajectory. Same parameters than in Fig. ( |2.4 ) 
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here, occurs when a system is periodically swept through an unstable state. The periodic 



amplification of fluctuations often leads to important consequences |49, 51]. 

A standard way of characterizing transient fluctuations of the stochastic process x{t) 
is by the time dependent moments {x'^{t)). This characterization gives the statistics of 
the random variable x at a fixed given time t. An alternative characterization is given by 
considering t as a function of x. One then looks for the statistics of the random variable t at 
which the process reaches for the first time a fixed given value x. The distribution of such 
times is the First Passage Time Distribution (FPTD). This alternative characterization 
emphasizes the role of the individual realizations of the process x{t). It is particularly well 
suited to answer questions related to time scales of evolution. For example, the lifetime 
of a given state can be defined as the Mean FPT (MFPT) to leave the vicinity of that 
state. The value of the associated variance of the PTD identifies whether that lifetime is 
a meaningful quantity. We will follow here the approach of first finding an approximation 
for the individual stochastic paths of the process, and then extracting the PT statistics 
from this description of the paths. In many cases the individual stochastic paths are 
easily approximated in some early regime of evolution, from which FPT statistics can be 
calculated. It is also often the case that statistical properties at a late stage of evolution 
can be calculated by some simple transformation of early time statistics. For example, we 
will see that there exists a linear relation between the random switch-on times of a laser 
calculated in a linear regime and the random heights of laser pulses which occur well in the 
nonlinear regime of approach to the steady state. 



3.1 Relaxational dynamics at a pitchfork bifurcation 

The normal form of the dynamical equation associated with a pitchfork bifurcation is 

dV 

dtx = ax + bx^ - cx^ + ^/e^(t) = -— + ^/e^(t), c>0 (3.1) 

ox 

where x{t) follows relaxational gradient dynamics in the potential V (see Section 4) and 
where we have added a noise term ^(f) with noise intensity e. The noise ^(i) is assumed to 
be Gaussian white noise of zero mean and correlation 



mat')) = 25{t-t'). (3.2) 

We are here interested in describing the relaxation triggered by noise from the state x = 0. 
In a supercritical bifurcation b < and we consider the relaxation for a > 0. In this case 
the quintic term in (3.1) is irrelevant and we can set c = 0. In a subcritical bifurcation 
b>0. When -f|^ < a < 0, x = is a metastable state and it becomes unstable for a > 0. 
Changing the control parameter a from negative to positive values there is a crossover in 
the relaxation mechanism from relaxation via activation to relaxation from an unstable 
state. The bifurcation diagram for the subcritical bifurcation is shown in Fig. (3.1). We 
will first consider the relaxation from x = in the supercritical case (a > 0,6 < 0, c = 0) 
and secondly the relaxation in the critical case a = of the subcritical bifurcation (6 > 0, 
c < 0). In the supercritical case the state x = is an unstable state and the relaxation 
process follows, initially, Gaussian statistics for x. However, the state x = has marginal 
stability in the case of a subcritical bifurcation for relaxation with a = 0, and in this 
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Figure 3.1: Bifurcation diagram for the subcritical bifurcation showing stationary solutions 
of X vs. the control parameter a. The form of the potential V for different values of a is also 
shown. Solid lines denote stable states. Dots and crosses stand, respectively, for unstable 
and metastable states, where cl = o-* = 

case Gaussian statistics, or equivalently, linear theory, does not hold at any time of the 
dynamical evolution. 

A usual approach to the statistical description of (3.1) features the Fokker Planck equa- 
tion (see Sect. 2.3) for the probability density P{x,t) of the process, (see Fig. (3.2)) 

dtPix,t) = d^{d^V)P + edlP = HP{x,t), (3.3) 

where H is the Fokker-Planck operator. Equations for the time dependent moments of the 
process x{t) are easily obtained from (3.3). The calculation of PT statistics is also standard 
lH, |5^ for a Markov process x(t): the probability that at time t, x(t) is still in an interval 
{xi,X2) given an initial condition xq at t = is known as the survival probability F{xo,t). 
It obeys the equation 

dtF{xo,t) = H+F{xo,t) (3.4) 

where is the adjoint of Hd. The latter operator differs from H in that it avoids the 
probability of reentering the interval (xi,X2) after having left it. The FPTD is given by 
f{xo,t) = —dtF{xQ,t), and as a consequence the MFPT, Tx^, obeys the Dynkin equation 

H+Tx, = -1. (3.5) 

The lifetime T of the state x = is given by the MFPT for to reach a value Rq starting 
at X = 0. From (3.5) we obtain 

T = l f""" dxie^(-^)A r e^(-^)Adx2. (3.6) 
e Jo Jo 
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This relation is generally valid. In the supercritical case, and for asymptotically small e, 



'(x,t) Jf) 



(b) 




(c) 





(d) 




Figure 3.2: Schematic evolution of the probability density P{x,t) for the decay from x = in 
a supercritical bifurcation: An initial distribution centered in x = broadens and develops 
two peaks around the deterministic steady state solutions x = ib(— Time evolution 
is (a) ^ (b) ^ (c) ^ (d). 



we obtain 



2a e 



2a 



— In e 

2a 



(3.7) 



where tp is the digamma function |1£]. The variance of the FPT distribution is given by 



(AT)^ 



4a2 



V^'(i). 



(3.8) 



We note that while in a deterministic treatment, x{t) starts to grow exponentially as 



^2at 



the lifetime of the state is not simply a~ because such lifetime is determined by 
fluctuations and it diverges logarithmically with e. On the other hand, (Ar)^ turns out 
to be independent of e in the limit of small noise intensity. The same type of calculation 



for the subcritical bifurcation |53, 54] gives a formula for T which interpolates between the 
Kramers result for relaxation from a metastable state, obtained for a < and ^1, and 
(3.7) for relaxation from an unstable state, obtained for a > and \k\ S> 1. The parameter 
k = a/(6e)^/2 measures the distance to the situation of marginality a = 0. 

As we stressed in the previous section much physical insight can be gained following the 
actual stochastic trajectories x{t) (Fig. (3.3)). We will now use this approach to describe 
the relaxation process and to reobtain the results quoted above for the distribution of 
passage times. A first intuitive idea of an individual trajectory in the supercritical case 
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Figure 3.3: A sample of individual trajectories for the decay process obtained from a nu- 
merical solution of (3.1) with x(0) = 0, a = 1, 6 = —1, c = 0, e = 10~*., and different 
realizations of the noise. The vertical bars indicate the difference between the most sep- 
arated trajectories at a fixed time. The evolution of an average trajectory {x^{t)) and a 
statistical spread 5{t) = (a:^(t)) — (|x(t)|)^ is also shown. The form of d{t) is associated with 
transient anomalous fluctuations. 

distinguishes two dynamical regimes in the relaxation process: There is first a linear regime 
of evolution in which the solution of (3.1) can be written (6 = c = 0) as 



where h{t) is a Gaussian process which plays the role of a stochastic time dependent initial 
condition which is exponentially amplified by the deterministic motion. There is a second 

time regime in which noise along the path is not important, and the trajectory follows 
nonlinear deterministic dynamics from a nonzero initial condition x(0). The deterministic 
solution of (3.1) (e = 0) is 



A theory for the relaxation process is based on this separation of dynamical regimes 
observed in each trajectory. In this theory the complete stochastic evolution is replaced by 
the nonlinear deterministic mapping of the initial Gaussian process h{t): The stochastic 
trajectory is approximated by replacing the initial condition x{Q) in (3.10) by h(t) so that 




(3.9) 



I- a 



x(0)e' 



(3.10) 



fx2(0)(e2«* - 1)]V2- 
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x{t), becomes the following functional of the process h{t): 

x(\h(t)],t) = . — ^^^^ -. (3.11) 

^ ^ [l-^/i2(t)(e2»*-l)]i/2 ^ ' 

Eq. (3.11) implies a dynamical scaling result for the process x{t), in the sense that x{t) is 
given by a time dependent nonlinear transformation of another stochastic process, namely 
the Gaussian process h(t). This result for x{t) is also often called the quasideterministic 
theory (QDT) |^] and has alternative presentations and extensions The fact that 

Eq. (3.11) gives an accurate representation of the individual trajectories of (3.1) (except 
for the small fluctuations around the final steady state Xy- = —a/b) justifies the approach. 

In the approximation (3.11) the PT are calculated as random times to reach a given 
value Rq within the linear stochastic regime. In fact, T sets the upper limit of validity of 
the linear approximation. The process h{t) has a second moment {h'^{t)) = |(1 — e"^'**) 
which saturates to a time independent value far from criticality, at ^ 1. In this limit the 
process h{t) can be replaced by a Gaussian random variable /i(oo) with {h(oo)) = and 
{h'^{oo)) = |. This permits to solve (3.9) for the time t* at which Rq is reached 

t* = —ln-^$^. (3.12) 

2a /i2(oo) ^ ^ 

This result gives the PT statistics as a transformation of the random variable h = h(oo). 
The statistical properties are, for example, completely determined by the generating func- 
tion 

W{X) = (e-^**) = / dhP{h)e-^''^^\ (3.13) 



where P{h) is the known Gaussian probability distribution of /i(oo). Moments of the PT 
distribution are calculated by differentiation with respect to A at A = 0. In this way one 
recovers (3.7) and (3.8). We note that this calculation assumes a separation of time scales 
between and T. When this separation fails there is no regime of linear evolution. The 
decay process is then dominated by fluctuations and nonlinearities (see an example in Sect. 
5.2). 

Time dependent moments of x{t) can also be calculated from the approximated trajectory 
(3.11) by averaging over the Gaussian probability distribution P{h,t) of the process h{t): 

(x"(t))= [ dhP{h,t)x''{[h{t)],t) (3.14) 



It turns out that in the limit e"* ^ 1 there is dynamical scaling in the sense that the 
moments only depend on time through their dependence in the parameter r = f e^"*. The 
result for {x'^{t)), given in terms of hypergeometric functions, can then be expanded in a 
power series in r |4^]. For example for the second order moment one obtains: 

oo 

{x\t)) = {x\^))Y.{-ir-\2n-l)\\T\ (3.15) 

n=l 

This result indicates that our approximation to the trajectory corresponds to a summation 
of the perturbative series in the noise intensity e which diverges term by term with time. 
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It also gives an interpretation of the MFPT as the time t for which the scahng parameter 
r ~ 1. For times of the order of the MFPT the amphfication of initial fluctuations gives 
rise to transient anomalous fluctuations of order e*^ as compared with the initial or final 
fluctuations of order e as shown in Fig. (3.3): At early and late times of the decay process 
different trajectories are close to each other at a fixed time, resulting in a small dispersion 
as measured by S{t). However, at intermediate times the trajectories are largely separated 
as a consequence of the amplification of initial fluctuations and 5{t) shows a characteristic 
peak associated with large transient fluctuations. 

The scaling theory discussed above can be put in a more systematic basis which can 
also be used in the subcritical case: In order to approximate the individual paths of the 
relaxation process deflned by (3.1) and x(0) = 0, we write x{t) as the ratio of two stochastic 
processes 

x{t) = z{t)/yy\t) (3.16) 
Then (3.1) (with c = 0) is equivalent to the set of equations 

dtz{t) = az{t) + ^^ey^'^{t)i{t) (3.17) 

dty{t) = -2bz\t), (3.18) 

with z{0) = x{0) = 0, y{0) = 1. Eqs. (3.17)-(3.18) can be solved iteratively from the initial 
conditions. In the zero-th order iteration 

x{t)r^z{t) = h{t)e''\ (3.19) 

where ^ 

h{t) = V^ [ e-^'^{s)ds, (3.20) 
Jo 

is a Gaussian stochastic process. In the first order iteration 

x(t) = ^"'^(^) (3 21) 

^ ' [l-26/o*e2«^/i2(s)ds]i/2- ^ ■ ^ 

In this approximation the decomposition (3.16) is interpreted as follows. The process z{t) 
coincides with the linear approximation (6 = 0,c = 0) to (3.1). The process y{t) introduces 
saturation effects killing the exponential growth of z(t). The scaling theory approximation 
(3.11) is recovered from this approach whenever at ^ 1 so that 

* e^'''h\s)ds ~ —h\t){e^^' - 1) (3.22) 
2a 

This indicates that the regime in which the scaling theory is valid is rapidly achieved far 
from the condition of criticality a ~ 0. 

We next consider the relaxation from x = with a = in a subcritical bifurcation 6 > 



|57|. As we did before we look for a decomposition of the form (3.16). Since we are only 
interested in the escape process we set c = in (3.1). Eq. (3.1) with a = c = is equivalent 
to the set (3.17)-(3.18) with a = 0. In the zero-th order iteration 



x{t) ~ z{t) = VeWit), W{t) = f\{s)ds. (3.23) 

Jo 
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The process x{t) coincides with the Wiener process giving diffusion in the locahy fiat po- 
tential V. In the first order iteration we find 



.it) = ^ -. (3.24) 

Two important differences between (3.24) and (3.21) should be stressed. The first one is 
that in the first stage of evolution given by (3.23) there is no escape from x = 0. The 
nonlinearities introduced by the process y{t) are essential for the escape process to occur. 
This means that, contrary to the case of (3.21), there is no regime of the escape process in 
which Gaussian statistics holds. The second difference is that x{t) given by (3.24) does not 
have a scaling form, since it is not a transformation of a single stochastic process. Indeed, 
x{t) in (3.24) depends on two non independent processes, namely W{t) and ^QW'^{s)ds. 
A naive approach to the problem would be to assume, by analogy with (3.11), a scaling 
form in which x{t) is given by the deterministic nonlinear mapping of the fiuctuating initial 
process ^/eW{t), i.e., x{t) given by the deterministic solution of (3.1) with x{Q) replaced 
by ^/(W{t). This would amount to take y{t) = 1 — 2ebtW'^{t). This scaling representation 
is qualitatively incorrect since it leads to a diverging MPT, T = oo. The accuracy of the 
representation (3.24) is shown in Fig. (3.4) where an individual trajectory given by (3.24) 
is compared with the corresponding one of the exact process given by (3.1). In this figure 
we observe that x{t) coincides initially with the Wiener process, it later departs from it and 
at a time rather sharply defined it departs from the vicinity of a; = 0. In fact, the strong 
nonlinearity implies that the solution of (3.1) with c = and x(0) ^ reaches |x| = oo in 
a finite time. It is then natural to identify the PT for the escape from x = as a random 
time t* for which |2;(t*)| = oo in (3.24), or equivalently y{t*) = 0. From (3.24) we find 

1 = 2be f W'^{s)ds. (3.25) 

Taking into account that W'^{es) = ew'^{s), (3.25) can be solved for t* as 

1 



2be9. 



1/2 



(3.26) 



with 



n= [^W^{s)ds. (3.27) 
Jo 

Eq. (3.26) gives the statistics of t* as a transformation of the statistics of another random 
variable $7. This scaling result for t* has the same basic contents than (3.12). In (3.12) 
the result appeared for at 1 as a consequence of Gaussian statistics while here the 
transformation (3.26) appears as an exceptional scaling at the critical point a = 0, and 
0, has non-Gaussian statistics. A discussion of the crossover between these two regimes is 



given in 1 53] 



The calculation of the statistics of t* from (3.26) requires the knowledge of the statistics 
of 0,. The latter is completely determined by the generating function G(A), for which an 
exact result is available []57| , 

G(A) = (e-^^) = (cosh(2Ai/2))-i/2_ (3,28) 
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Figure 3.4: Stochastic paths in the decay from x = for a = in the subcritical bifurcation 
with b = c = 1 and e = 10~^. Different paths correspond to different approximations for 
the same reahzation of the noise. The sohd Une corresponds to a simulation of the exact 
process given by (3.1). Dotdashed hne corresponds to a Wiener process. Dots correspond 
to the scaUng approach similar to (3.11) (see text) and the dashed line corresponds to the 
approximation (3.24). 



The moments of the PTD are obtained from (3.26) in terms of G{X) as 



(f = " / dX A"-iG(A). (3.29) 

r(n) Jo 

The transient moments {x"'{t)) can also be obtained in terms of the statistics of the PT 
|57|: Given the strong nonlinearity of the escape process a good approximation to calculate 



ensemble averages is to represent x{t) by 

x'^(t) =x'^{oo)e{t-t*), (3.30) 

where x{oo) is the final stable state (local minima of the potential V) and 9{t — t*) is 
the Heaviside step function. The individual path is approximated by a jump from x = 
to a; = x(oo) at a random time t* . The transient moments are then easily calculated as 
averages over the distribution of the times t* . 

The methodology described here to characterize relaxation processes by focusing on the 
statistics of the passage times can be generalized to a variety of different situations which 
include colored noise |58], time dependent control parameter a{t) sweeping through the 



instability at a finite rate |59, p0[, competition between noise driven decay and decay 



induced by a weak external signal |61], role of multiplicative noise in transient fluctuations 
ll62|| , saddle-node bifurcation etc. We have limited here ourselves to transient dynamics 
in stochastic processes defined by a SDE of the Langevin type, but the PT characterization 
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is also useful, for example, in the framework of Master Equations describing transport in 
disordered systems |54, ^ 



3.2 Statistics of laser switch-on 

Among the applications to different physical systems of the method described above, the 
analysis of the stochastic events in the switch-on of a laser is particularly interesting for 
several reasons. Historically (see the review in [^]) the concept of transient anomalous 
fluctuations already appeared in this context in the late 60's, and from there on a detailed 
experimental characterization exists for different cases and situations. The idea of noise 
amplification is also very clear here, and for example the term "statistical microscope" 
IgsII has been used to describe the idea of probing the initial spontaneous emission noise 
through the switch-on amplification process. Prom the methodological point of view the 
idea of following the individual stochastic trajectory has in this context a clear physical 
reality. While in more traditional equilibrium statistical physics it is often argued that 
measurements correspond to ensemble averages, every laser switch-on is an experimental 
measured event which corresponds to an individual stochastic trajectory. Also here the 
PT description has been widely used experimentally |^6[. Pinally, from an applied point 
of view, the variance of the PT distribution for the switch-on of a semiconductor laser sets 
limitations in the maximum transmission rate in high speed optical communication systems 
IgT] , |68|| . We will limit our quantitative discussion here to semiconductor lasers. A review 
of other situations is given in |48|. 

The dynamics of a single mode semiconductor laser can be described by the following 
equations for the slowly varying complex amplitude of the electric field E and the carrier 
number |9|, 0, |7l|, 



1 -\- f fy 

dtE=——[G--f]E+^CEit), (3.31) 



dtN = C- - G\E\^ - V4pN[E*^E{t) + ECE{t)]. (3.32) 

Terms proportional to the gain coefficient G account for stimulated emission and, to a first 
approximation, G = g{N — Nq), where 5 is a gain parameter and A^o the value of carrier 
number at transparency. The control parameter of the system is here the injection current 
G. Spontaneous emission is modeled by a complex Gaussian white noise S,Eit) of zero mean 
and correlation 

{CE{t)CE{t')) = 26{t-t'). (3.33) 

Noise is here multiplicative and its intensity is measured by the spontaneous emission rate 
p. The noise term in the equation for A^ can be usually neglected 0. The variables E and 
A have decay rates 7 and 7e such that 7 ~ 10^ 7e. This implies that the slow variable is A, 
and that when the laser is switched-on the approach to the final state is not monotonous. 
Rather, a light pulse is emitted followed by damped oscillations called relaxation oscillations. 
This also occurs for GO2 and solid state lasers, but there are other lasers such as He-Ne 

*When the noise term in (3.32) is neglected, the Ito and Stratonovich interpretations are equivalent for 
(3.31) 
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("class Alasers") in which E is the slow variable, the approach is monotonous and a simpler 
description with a single equation for E holds. Such description is similar to the one in Sect. 
3.1 I^Sj. Finally, the a-parameter, or linewidth enhancement factor, gives the coupling of 
phase 4> and laser intensity I {E = I^^'^e'^'^). 



dtI=iG-j)I + ipN + ./IpNI^iit), (3.34) 



9*<^=f(G-7) + y^e<^(0, (3.35) 

where and are the intensity and phase noise components of (,e- These equations are 
conventionally written as Ito SDE. The instantaneous frequency of the laser is the time 
derivative of the phase. Transient frequency fluctuations for class A lasers were discussed 
in@. 

To switch-on the laser, the injection current C is switched, at a time t = 0, from below to 
above its threshold value Cth = ( ^ +-^o)7e. Below threshold the laser is off with / fluctuating 
around 1 = 0. The steady state value above threshold is given by I = ^{C — Cth) + and 

N = Nth = • An example of the early time dynamics in the approach to this steady state 
is shown in Fig. (3.5) where two switch-on events corresponding to two different stochastic 
trajectories are shown. We observe that the laser switches-on when N reaches a maximum 
value. This happens at a random time which gives the delay from the time at which C is 
switched. The dispersion in the delay times is known as "jitter". We will characterize the 
switch-on statistics by the PT distribution for / to reach a given reference value. We also 
observe a statistical spread in the height Im of the light pulses (maximum output power) 
for different switch-on events, being larger for longer delay times. The frequency of the 
laser shows huge fluctuations while it is off, and it drifts from a maximum to a minimum 
frequency during the emission of the light pulse due to the phase-amplitude coupling caused 
by the a-parameter. This excursion in frequency is known as "frequency chirp", and again 
there is a statistical distribution of chirp for different switch-on events. Relevant questions 
here are the calculation of switch-on times, maximum laser intensity and chirp statistics, 
as well as the relation among them. 

The calculation of switch-on time statistics follows the same basic ideas than in Sect. 3.1 
[70, 72]: We consider the early linearized regime formally obtained setting G = in the 



equation for A^. When the solution of this equation is substituted in (3.31) we still have a 
linear equation for E but with time dependent coefficients, so that, 

E{t) = e^^'^/^h{t), (3.36) 



A{t) = (1 + ia) j'_ dt'[g{N{t') - N^) - 7], (3.37) 

h{t) = J_' dt'^(3N{t')^{t')e-^^''^/^, (3.38) 

where we have chosen a reference time t, N(t) = Nth- We note that for ordinary parameter 
values N{t) grows linearly with t around t (as it is seen in Fig. (3.5)), so that A{t) ~ {t — t)"^. 
The important fact is again that the complex Gaussian process h{t) saturates for times of 
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Figure 3.5: Time evolution of (a) laser intensity / = \E\'^ , (b) carrier number and (c) 
frequency / = {2ii)~^^ as obtained from the numerical integration of Eqs. (3.31)-(3.32) 
for typical parameter values and for two different switch-on transients: solid and dashed 
line. Times are given in picoseconds. 



interest to a complex Gaussian random variable h(oo) of known statistics. We can then 
solve (3.36) for the time t* at which the reference value Ir = \Er\'^ of / is reached 

The statistics of t* are now easily calculated as a transformation of that of h{oo). In 
particular we obtain (AT) <^ T, with 



{n=t + J J l-^], r = ln:^ (3.40) 



g{C - Ctn) ' 2r 
2Tg{C - Cth) 



(^AT)^ = (3.41) 



where b is proportional to /?. 

We next address the calculation of the statistics of the maximum laser intensity Im 
|)70| , [T^ . Fig (3.6) shows a superposition of different switch-on events, obtained from a 
numerical integration of (3.31)-(3.32), in which the statistical dispersion in the values of Im 
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Figure 3.6: Superposition of 50 switch-on events (50 individual stochastic trajectories) for 
the laser intensity as obtained from the numerical integration of Eqs. (3.31)- (3.32). The 
laser power P is proportional to the laser intensity I. 



is evidentiated. The laser intensity peak value is reached well in the nonlinear regime and 
we adopt, as in Sect 3.1, the strategy of mapping the initial stochastic linear regime with 
the nonlinear deterministic dynamics. For each switch-on event, the deterministic equations 
should be solved with initial conditions at the time t* which identifies the upper limit of the 
linear regime: E = E{t*) = Er and N = N{t*). The value of N{t*) is calculated hnearly 
and takes random values for different switch-on events. Eliminating the parameter t the 
solution can be written in an implicit form as 

E = E{N,Er,N{t*)). (3.42) 

Any arbitrary function Y of the dynamic variables in which one might be interested can 
then be written as 

Y{t) = Y{N{t),E{t)) = Y{t, N{t*), Er) = Y{N, Er, N{t*)). (3.43) 

To proceed further one generally needs an explicit form for the solution of the equations. 
However, important results can be already obtained if one is interested in extreme values 
of Y. The variable Y takes an extreme value Y^ at a time such that dtY{t) = 0. This 
gives an additional implicit relation N{te) = N(.{Er, N(t*)) so that, 

Y{te) = Ye{Er,N{t*)). (3.44) 

This equation gives the desired mapping of the statistical properties of the switch-on times 

t* into those of an extreme value of a dynamical variable: For each switch-on event there is 
a value of t* and a corresponding value Yg(t*). Taking advantage of the fact that (AT) <C T 
we can expand Ye{t*) around t* = T: 



Y(n-Y(T) + ^^^ 



{t* -T) + 0{ATf. (3.45) 

T 
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Applying this relation to the maximum output power, = Im, it predicts a straight 
line for a plot of the values of Im vs. switch-on times obtained from many independent 
switch-on experiments. Such linear dependence can be already seen in Fig. (3.6). This 
prediction has been corroborated experimentally for CO2 0, |7^, semiconductor and 
solid state lasers [^]. The slope of this linear relation can be estimated in different cases. 
For semiconductor lasers |7Cl| j it is easily seen from the equations above, that the slope is 
positive and proportional to the distance to threshold C — Cth- However in other cases 
there is a predominant contribution of the initial fluctuations of N{t = 0) and the slope 
becomes negative |74]. 

Turning now to transient frequency fluctuations, a first delicate issue is the proper defi- 
nition of an instantaneous frequency, since dt(j) in (3.35) contains a white noise noise contri- 
bution, and therefore, in principle, an infinite variance. A meaningful operational definition 
of the instantaneous frequency can be given in terms of a short time interval average of dt4), 
or equivalently in terms of a time resolved field Fourier spectrum |71, 77|. From (3.35) we 
also see that there are two sources of frequency fiuctuations. One is the noise along the path 
given by the phase noise ^<^, and the second is the stochastic nature of the gain term G for 
different switch-on events, due to the random nature of t* . The first dominates steady state 
properties such as laser linewidth, but the second dominates transient properties such as 
the fluctuations in the chirp range. Since the chirp range is defined by a difference between 
two extreme values of the frequency the general framework discussed above applies. In fact 
one finds that the maximum and minimum instantaneous frequencies during the transient 
are linear functions of the switch-on time t* (T^ . Chirp noise statistics is again calculated 
by a simple transformation of the switch-on time statistics. 

We finally mention that the operation of a semiconductor laser in a system of optical 
communications involves a repetition of switch-on events (large signal modulation) to trans- 
mit pulses at a given bit-rate. Communication errors are likely to happen when AT is large 
or when memory effects are present. The effect of memory occurs when a given switch-on 
event is not independent of the previous one in very fast modulation. The determination of 
" memory- free" regions in parameter space for a given modulation rate involves, as a crucial 
time scale of the problem, the mean switch-on time T |7|. Other errors are due to pulse 
dispersion in the optical fibre because of large chirp noise [79|. 



4 Noise in Dynamical Systems 

4.1 Classification of Dynamical Systems: Potentials and Lyapunov Func- 
tions 

In this section, we will study the effect of noise in the long time behavior of dynamical 
systems. First of all, though, we will briefiy review the basic concepts about the stability 
of dynamical systems ||80t| . 

Let X = (xi, . . . , xat) be a set of real dynamical variables. They follow the evolution 
equations: 

dx ' 

-^ = fi{xi,...,XN), i = l,...,N (4.1) 
One is usually interested in finding the fixed points x of the dynamical system as those 
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Figure 4.1: Illustration of the concepts of stable and asymptotically stable fixed points. 



points which do not evolve in time, i.e. those satisfying: 

= 0, i = l,...,N (4.2) 



dxi 



or, equivalently, 

/,(xi,...,X7v) = 0, i = l,...,N (4.3) 

Knowing the fixed points is as important as knowing their stability. Loosely speaking, a 
fixed point x is stable if any initial condition x(0) sufficiently close to x remains close to x 
as time evolves. The fixed point is said to be asymptotically stable if any initial condition 
x(0) sufficiently close to x tends to x as time evolves. Fig. (4J.) illustrates the concepts of 
stability and asymptotic stability. 

To determine the stability of a fixed point, Lyapunov's stability theorem is of great help. 
The theorem states that if we can find a Lyapunov function, L{x) = L{xi, . . . ,xn) which 
is such that it has a local minimum at x = x and monotonically decreases with time, i.e.: 

(a) L{x) > L{x) for X 7^ X 

dL ^dLdx, (4.4) 
(6) — = > — < for X 7^ X 

^ ' dt dxi dt - ^ 

then X is a stable fixed point. If condition (b) is replaced by 

(6') ^ < for X / X (4.5) 
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then X is an asymptotically stable fixed point. L{x) is also called a Lyapunov potential or, 
simply, the potential. 

The description of the behaviour of a dynamical system is greatly simplified if we know 
of the existence of a Lyapunov potential: the system evolves towards the minima of the 
Lyapunov and once there it stays nearby when small perturbations act upon the system. 
The Lyapunov potential not only tells us about local stability, but also gives information 
about global stability. If two fixed points x^-^^ and x*-^-* are such that L{x^'^^) < L{x^^^) then 
we can infer that x^^) is the stable fixed point and that x^^^ is a metastable fixed point. 
This means that a sufficiently strong perturbation might take the system out of x^^^ to 
x^'^\ Therefore, to understand the asymptotic dynamics it is of great importance to find 
out whether a Lyapunov potential can be found for the dynamical system under study. It 
is not an easy task to (i) prove that the Lyapunov potential does exist, (ii) find it. The 
systems for which a Lyapunov function exists are called in the literature potential systems 
and the name non-potential systems is therefore reserved for those systems which do not 
have a Lyapunov function (notice that this is different from the fact that there is no known 
Lyapunov potential). According to these definitions, and following the ideas of Graham on 



nonequilibrium potential |81, 32, 84], the following classification of dynamical systems 

can be established |8^ : 

1.- Relaxational Gradient Flow 

2 - Relaxational Non-Gradient Flow 

3. - Non-Relaxational Potential Flow 

4. - Non-Potential Flow 

The first three categories are called potential flows. We now describe each of them in some 
detail. 

(1) Relaxational Gradient Flow. 

Those are systems for which there exists a function (the "potential") V{x) = V{xi, . . . ,xm) 
in terms of which the dynamics is written as: 

The fixed points of this dynamical system are the extrema of the potential V{x). The 
trajectories lead to any of the minimum values of V{x) following the lines of maximum 



slope (steepest descent) as indicated schematically in Fig. ( |4.2| ). In this case the potential 
y is a Lyapunov function. The proof of this statement is very simple indeed (we need the 
additional condition that V{x) is bounded from below): 

^=f ^^ = -f f^V<0 (4.7) 
dt ^ dxi dt ^ \dxiJ - ^ ' 

An example of this class of systems is the real Ginzburg-Landau equation (see later), 
also known in the theory of critical dynamics as model A in the classification of (see 
Section 6). 

(2) Relaxational non Gradient Flow. 

There exists again a potential function, V{x), but now the dynamics is not governed just 
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Figure 4.2: Schematic trajectories for the case of relaxational gradient flow in two cases: 
a Lyapunov function with a single minimum (left picture) and a Lyapunov function with 
two minima (right picture). The system evolves towards one of the two minima, according 
to initial conditions, following the lines of maximum slope. 



by V, but it is given by: 



dxi ^ dV 



dt dxj 



where Sij{x) is a real, symmetric, positive definite matrix. The fixed points of the dynamics 
are still given by the extrema of V and the trajectories lead to the minima of V but not 
necessarily through the lines of maximum slope. In this sense, we can say that the transient 
dynamics is not governed just by the potential. Fig. (^^) shows schematically the situation. 
A typical example of this dynamics is the so-called Cahn-Hilliard |p8| model for spinodal 
decomposition or model B in the context of critical dynamics [^]. In this model the 
dynamical variables represent a concentration field that follows a conservation law: 

, N 

= (4.9) 

1=1 

This constrain is satisfied by choosing 5 = — V^, the lattice Laplacian as defined, for 
instance, in ( 2.134] ). 



In this case, the potential V is still a Lyapunov functional: 



dt ^ dxi dt ^ dxi dxi - ^ ^ 

1=1 * jj=i * 3 
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Figure 4.3: Schematic trajectories for the case of relaxational non-gradient flow. The 
system evolves towards the minimum of the potential, but not necessarily through the lines 
of maximum slope. 

since, by definition, S* is a positive definite matrix. 
(3) Non Relaxational Potential flow. 

These systems are characterized because the asymptotic behavior is not determined simply 
by the minima of the potential, but there is a residual dynamics once those minima have 
been reached. A first category within this class is given by: 

Here, Dij{x) is an arbitrary matrix that can be split into symmetric and antisymmetric 
parts: 

and we demand that S is positive definite matrix. Again, the fixed point of the dynamics 
are determined by the extrema of the potential V and y is a Lyapunov potential: 

— - - V 5 - V A < f4 13) 

dt ^■^ dxi dxi dxi dxi ~ 

1,3=1 ' J «j=i 

The second sum of this expression is zero due to the antisymmetry of matrix A. We are 
then led to the situation sketched in Fig. (^^) : the system evolves towards the minima of V 
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Figure 4.4: Schematic trajectories for the case of non-relaxational potential flow. The 
system evolves towards the (degenerate) minima of the potential V and, once it has reached 
there, there is a residual movement on an iso-potential region. 

but once it reaches there, there is a residual movement produced by the antisymmetric part 
of the dynamical coefficients that can make the system flow without cost of the Lyapunov 
function. The dynamics of the nematodynamics equations commonly used to describe the 
dynamics of liquid crystals in the nematic phase belongs to this class, as shown in [^9| . 

The second, more general, category of non-relaxational potential flow will be one in 
which the dynamical equations can be split into two parts, namely: 



Here S{x) is a symmetric, positive definite matrix and Vi{x) represents the residual dynamics 
after the relaxational part has acted. Since we want F to be a Lyapunov potential the 
residual dynamics must not contribute to its decay: 



dxi 



N 



dV 




(4.14) 




(4.15) 



or, since the first term of the r.h.s. is always negative, a sufficient condition is: 




(4.16) 
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This is the so-called orthogonality condition. By using (4.14) can be written as: 

1=1 \ j=i 

or, in compact vector notation: 

(/ + svv) • vy = 




(4.17) 



(4.18) 



It can be shown that this equation has the structure of a Hamilton-Jacobi equation and in 
general it is very difficult to solve to find the potential V, although important results have 



been obtained in the literature |81 



(4) Non Potential flow. 

Finally, a non-potential flow is one for which the splitting (|4.14| ), satisfying ( 4.15| ), admits 
only the trivial solution V = 0, Vi = fi. Notice that for a system to be classified as non- 
potential, we have to prove the non-existence of non-trivial solutions of ( [4.17| ) which, of 
course, it is not the same that not being able to find non-trivial solutions of (|4.1'i 



Let us consider in detail an example of the previous classification scheme: the complex 
Ginzburg-Landau equation (CGLE) [ 



91|, |2|. 



Let A{f, t) be a complex field satisfying the equation of motion: 

P\A\^A 



dA , ^2 . 



(4.19) 



= /^i? + V/j P = Pr + if^i and 7 = 7i? + *7/ are complex numbers. This is the so-called 



complex Ginzburg-Landau equation Q. It appears in a great variety of situations and 
it is a generic equation for the slow varying amplitude of the most unstable model near a 
Hopf bifurcation in an extended system (see Sect. 5.3). This might be put in what is called 
in the literature a "potential" form: 

dA _ _5V_ 
'dt ~ 



SA* 



(4.20) 



where the "potential" V is given by: 
V{A,A*} -- 



dr 



-/.|^|2 + 7|VA|2 + ||>1|4 



(4.21) 



This, being of course true, is somewhat misleading because the function V does not have 
most of the properties associated with a potential function as defined in the previous sec- 
tions. In particular it can not be a Lyapunov function because it generally takes complex 
values. 

We pose now the problem in a different form. Considering the evolution of the real and 
imaginary components of A{f,t) = X[f,t) + iY{f,t), is it possible to find a real function 
V{X, Y) such that the flow is relaxational gradient? i.e. 



(4.22) 



dX 


5V 


~dt 


~6X 


dV 


6V 


'dt 
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This would imply: 

A /rlY\ A /Piv\ 

(4.23) 











5Y 


if) 


5X 


(f) 


* 


7 


* 

= 7 


13 = p* 



which leads to 

(4.24) 

i.e. a complex Ginzburg-Landau equation with real coefficients. 

In general, though, when the coefficients are complex there is no obvious Lyapunov 
function. However, there is an interesting case in which the dynamics can be put in the 
category (3) described above: there exists a split of the dynamics in terms of a potential 
term plus an extra term which, being orthogonal to the first, does not contribute to the 
variation of the Lyapunov potential. The easiest way of showing this is noticing that if we 
write V{A,A*) = Vr{X,Y) + iVi{X,Y) the variation of the real part is: 



dt 2 J 



(&Vr\'^ , (&Vi 
\5X ) ^ \5Y 



+ 



(4.25) 



2(7//9/? - /3/7r) / dr{X^ + Y^)V{YVX - XVY) 
It is clear now that if the following condition is fulfilled 

7//3i? - (illR = (4.26) 
Vr acts as a Lyapunov. Introducing V = V/2 the dynamical system can be written as: 



dx 
dt 
ax 
dt 




(4.27) 



and the orthogonality condition is: 




= (4.28) 

which is satisfied by condition ( [4. 261 ) . In this case a change of variables of the form A{f, t) 
A{f, t)e*°* for a convenient value of a leads to 

(-) 

The 1 in the r.h.s. is the real CGLE and the i in the r.h.s. gives rise to a form of the 
Nonlinear Schrodinger equation. 

We finally mention that although for general complex coefficients fi, (3 and 7 there is 
no simple form for a Lyapunov functional, this does not mean that there is no solution 
for the orthogonality condition ( [4.17| ). In fact, Graham and collaborators [91, ^ have 



found an approximate Lyapunov function for the CGLE by a gradient expansion solution 
of the orthogonality condition ( 4.17| ). A numerical test of its Lyapunov properties and a 



discussion of such Lyapunov function for the CGLE is given in 
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4.2 Potentials and Stationary Distributions 

We will see now which is the effect of noise on the dynamics of the above systems. To this 
end, we consider the general dynamics with the addition of a general noise term (to be 
considered in the Ito sense): 

N 

Xi = fi{x) + Y,gijij{t) (4.30) 
i=i 

where gij{x) are arbitrary functions and ^j{t) are white noise variables: random Gaussian 
variables of zero mean and correlations: 



(ei(i)ei(i')) =2e<5i,'^(t-t') 



(4.31) 



e is the noise intensity, which is included here explicitly. Due to the presence of the noise 
terms, it is not adequate to talk about fixed points of the dynamics, but consider instead 
the maxima of the probability distribution function P(x,t). This satisfies the multivariate 
Fokker-Planck equation m: 



dP{x,t) 
di 



N 



N 



d 



7 = 1 ^ 



where the matrix G is 



(4.32) 



(4.33) 



For the case of relaxational (gradient or non-gradient) flows, an explicit solution in the 
asymptotic regime of the Fokker-Planck equation can be found if the relation G = S \s 
satisfied and 5 is a constant matrix. In this case one can show easily that the stationary 
solution is given by: 

V{x)' 



Pst{x) = iVexp 



(4.34) 



In some sense, we can say that the stationary distribution is still governed by the minima 
of the potential, since the maxima of this stationary distribution will be centered around 
these minima, see Fig. (|4.5| ). In Statistical Mechanics we encounter a large variety of 
dynamical systems that can be described by the stationary distribution ( f4.34| ) []87| . The 
potential V is the free energy F of the system and the noise terms have a thermal origin 
and, consequently, the noise intensity is nothing but the absolute temperature, e = fesT. 
The equality G = S* is known as the fluctuation-dissipation relation and it is the one 
that ensures that the asymptotic (equilibrium) distribution of a dynamical model deflned 
by eqs. (4.30) and ( [4.31| ) is the Gibbs distribution exp[— F/ZcbT"]- Sometimes, noise has 
to be added to the equations of motion in order to reach thermal equilibrium. This is, for 



instance. Cook's modiflcation [93| to the Cahn-Hilliard model for spinodal decomposition 



mentioned above |46|. 

The situation for non~relaxational flows is somewhat more complicated. General condi- 



tions for the validity of ( 4.34 ) and their relation with propoerties of detailed balance have 
been discussed in [p4|. We explain now some particular results. Let us consider a general 
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Figure 4.5: Potential and stationary probability density function 



non-relaxational potential flow with the addition of noise terms: 



N 



(4.35) 



with 



(4.36) 



The noise variables ^j{t) are Gaussian distributed with zero mean and correlations given 



by ( 4.31 ). It can be easily shown that: if (i) 5 is a constant matrix, (ii) the fluctuation- 
dissipation relation G = S holds, (iii) the residual term vi satisfies the orthogonality con- 
dition, eq. ( 4.16| ), and (iv) Vi is divergence free: 



^ dxj 







(4.37) 



then the stationary distribution is still given by ( 4.3^ ). In this case Vi can be easily in- 
terpreted as a non dissipative contribution to the dynamics since it does not contribute to 
the time variation of V and ( 4.37D implies that it conserves phase space volume. These 
conditions are satisfied, for example, if the non-dissipative contribution is of the form 
Vi = — J2j ^ijdxjV with A an antisymmetric matrix. 

An important result, due to Graham pi], 83, pi] states quite generally that: For a 



non-relaxational potential flow (i.e. the orthogonality condition ( [4.17 ) is satisfied) and 
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the fluctuation-dissipation relation S = G holds, then the stationary distribution is still 
governed by ( 4.34 ) in the small noise limit. More precisely, we have 



V{x) = lim -e In P,t(x) (4.38) 

In those case the potential V{x) is given also the name of "Graham potential", and the 
effect of the noise terms on the asymptotic dynamics is to introduce fluctuations (governed 
by the Gibbs distribution) around the remanent dynamics which occurs in the atractors 
identified by the minina of V. 



4.3 The Kiippers-Lortz instability 

When a system displaying the Rayleigh-Benard instability is put under rotation, the 
action of the Coriolis forces induces what is called the Kiippers-Lortz (KL) instability [p6| . 
This is relevant to many geophysical and astrophysical fluids which are subjected to, e.g., 
planet rotation. The KL instability can be described as follows: for an angular rotation 
speed, ri, greater than some critical value Oc, the convective rolls, in the reference frame 
that rotates with the system, alternatively change orientations between three preferred 
directions. 

Based on the fact that, in a first approximation, only three directions are relevant to 
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this problem, Busse and Heikes |97] have proposed a dynamical model to study the KL 
instability. In this model the amplitudes of the three rotating modes, Ai, A2, A3, follow 
the evolution equations: 



(4.39) 



// is a parameter related mainly to the temperature gradient and 6 is related to the rotation 
speed in such a way that = implies 6 = 0. We will consider only 17 > 0, or 5 > 0; the 
case < (5 < 0) follows by a simple change of the coordinate system. Finally, although 
Ai are complex numbers, their phase effectively disappears from the previous equations and 
they can be considered real variables. A similar set of equation has been proposed to study 
population competition dynamics. If we a have a single biological species, the Verhulst or 
logistic model assumes that its population N{t) satisfies the evolution equation: 

dN 

— = rN{l - aN) (4.40) 

where r is the reproductive growth rate and a is a coefficient denoting competition amongst 
the members of the species (for food, for example). If three species are competing together, 
it is natural to modify these equations to a Gause-Lotka-Volterra form [98|: 

iVi = riA^i (1 - A^i - aN2 - /SN^) 

N2 = r2N2il- N2-aN3- (3Ni) (4.41) 
iVs = rsNs{l-Ns-aN,-pN2) 

which, after suitable rescaling, are the same that the Busse-Heikes equations for the modulus 
square Oj = \Ai\'^ of the amplitudes. 
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H 



Figure 4.6: Stability regions for the Busse-Heikes dynamical system ( [4.39| ). The region H 
is where the hexagon solution (three equal amplitudes) is stable. In the R region, the three 
roll solutions are stable, and in region KL there are no stable fixed points. 

The fixed point solutions of the Busse-Heikes equations are either of two types: (a) roll 
(R) solutions and (b) hexagon (H) solutions. 

(a) Roll solutions. There are three of these solutions, each characterized by a unique 
non-vanishing amplitude, for instance: (^1,^2,^3) = (1,0,0) is a roll solution with rolls 
parallel to the ei direction, and so on. 

(b) Hexagon solutions. These are stationary solutions in which the three amplitudes are 
different from 0, namely Af = A2 = A"^ = 3^^- 

The stability of the roll and the hexagon solutions can be studied by means of a linear 



stability analysis. The result is summarized in Fig. (4.6). For fi < the only stable situation 
is the hexagon solution. For ^ > 0, but fi > 5, the roll solution is stable and when fj, < 5 
however, the roll solution is unstable. This is the KL instability that occurs at < /x = (5 
and can be described as follows: consider the roll solution {Ai, A2, A3) = (1, 0, 0), the mode 
A2 starts growing and Ai decreasing in order to reach the roll solution (0, 1,0). However, 
this new roll solution is also unstable, and before it can be reached, the dynamical system 
starts evolving towards the roll solution (0,0,1), which is unstable and evolves towards 
the solution (1,0,0) which is unstable ... This is the KL instability that shows up in the 
rotation of the convective rolls. 

(1, 0, 0) ^ (0, 1, 0) ^ (0, 0, 1) ^ (1, 0, 0) ^ (0, 1,0)... (4.42) 

In brief: we have in the KL region 3 unstable fixed points, each one of them evolves in the 
next one. 
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Further insight into the dynamical system can be obtained by rewriting the BH equations 
showing up explicitly that the terms proportional to S do not have a relaxational origin. If 
we introduce the potential function V{Ai,A2,A3): 

V{A,,A2,As) = =^{Al+Al+Al) + -^{At+At+Al) + ^{AlAl+AlAl + AlAl) (4.43) 
we can write: 

where 

VI = 6Ai{-Al + Al) 

V2 = 6A2i-Al + Al) (4.45) 

V3 = 6A3{-Al + Al) 

If (5 = the system is relaxational gradient and the corresponding stability regions can 
be obtained also by looking at the minima of the potential V. For the roll solution, the 
potential takes the value: 

Vr = ^ Roll solution (4.46) 

whereas for the hexagon solution: 

-3 

Vh = 77 r Hexagon solution (4.47) 

4(3 + 2fi) 

According to this potential function, and for 5 = 0, the rolls are more stable than the 
hexagons whenever Vr < Vh, i-e. /i > 0. Unfortunately, this simple criterion does not have 
an equivalent in the non relaxational case, 5^0. 

The potential V acts as a Lyapunov function in the potential case, 6 = 0. This is 
obvious from the form ( 4.44| ) of the dynamical equations. For the case 6 ^ we write down 



the orthogonality condition ( [4.16| ) and it is satisfied by the same potential function ( 4.43| ) 
provided the following condition is satisfied: 

dfi = (4.48) 

The case (5 = is the relaxational gradient case. The case = is very interesting because 
the movement appears to be non-relaxational potential with a known potential. We will 
explore later this interpretation. Now we solve "essentially" the dynamical equations for 
fi = 0. By "essentially" we mean that we can establish the asymptotic period of movement, 
after a transient time of order 1. First we write down the general equations for the variables 
tti (time has been rescaled by a factor of 2): 

(ii = ai{l - cii - 1 + fi{a2 + as) - 6{a2 - a^)) 

(12 = a2(l - ai - 1 + /i(a2 + as) - (5(a3 - ai)) (4.49) 

da = 03(1 - oi - 1 + //(a2 + as) - 5(ai - 02)) 

If we introduce now the variable x{t) = ai + 02 + 03, it is straightforward to show that it 
satisfies the evolution equation: 

X = x{l - x) - 2fiy (4.50) 
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where: 

y{t) = aia2 + 0203 + a^ai (4.51) 
In the case fi = the equation for x{t) is a closed equation whose solution is 

x{t) = -J- \— (4.52) 

(here = x(i = 0)). From this Gxprcssion it turns out that liiii^_i.cxD 

x{t) = 1 independently 

of the initial condition. In practice, and due to the exponential decay of the above expres- 
sion, after a transient time of order 1, x{t) already takes its asymptotic value x{t) = 1. 
Therefore, we can substitute ai(t), say, by 1 — 02 (i) — o,^{t) to obtain evolution equations 
for 02 (i) and a^{t). In this way, the original 3- variable problem is reduced to a 2-variable 
one: 

0,2 = -(502(1-02-203) f4 53) 

03 = -(5o3(l - 202 - 03) 

These equations have a Hamiltonian form: 



n — j^9n 

«2 - 075777 

with the Hamiltonian: 



Tan (4.54) 



7^ = 0203(1-02-03) (4.55) 

As a consequence of this Hamiltonian structure for the dynamics valid for fi = 0, it turns 
out that the "energy" H = Ti.{t) is a constant of movement after the transient regime of 
time (of order 1). After this transient, the Hamiltonian description is valid. The movement 
in the line Tl{t) = H, see Fig. ( |4.7| ), is periodic in the variables 02, 03 (and hence in the oi 
variable also). 

The value of H depends only on initial conditions. To check this let us use the following 
equivalent (in the asymptotic time limit) definition of TC: 

TL = 010203 (4.56) 

From this definition, and for arbitrary value of /i, we obtain the following exact equation 
for H: 

HI-/ 

n~^— = 3- {3 + 2fi)x (4.57) 

(one can reduce the original Busse-Heikes model to dynamical equations for variables 
x, y, Ti but the equation for ij turns out to be too complicated, see |Q). If we sub- 
stitute the solution ( [4.52| ) for x{t) valid in the case = we obtain: 

7^(t) =Ho [(l-xo)e~* + xo]^^ (4.58) 
with TLq = 7i{t = 0) = 01(0)02(0)03(0). The asymptotic value for 7i is 
H - lir^ -Nff) - ->7 ^-3 _ ai(0)o2(0)o3(0) 



(Oi(0)+O2(0)+O3(0))2 
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Figure 4.7: Asymptotic dynamics for = 0. Here we plot the level curves 0203(1 — 02 — 
as) = H for the Hamiltonian Ti. defined in ( [4.55| ). The allowed values for H range between 
H = 1/27 (most inner level curves) and H = external level curves. 



Again, this asymptotic value is reached after a transient time of order 1. 

The behavior of the dynamical system in the case /i = can now be fully described 
and interpreted: after a transient time (or order 1) the three variables oi, 02, 03 vary 
periodically in time such that 01 + 02 + 03 = 1. When oi decreases, 02 increases, etc. 
This is characteristic of the Kuppers-Lorz instability. The motion is periodic because it is 
a Hamiltonian orbit with a fixed energy. The orbits can be plotted as closed trajectories 
in the 01 + 02 + 03 = 1 plane. The exact shape of the trajectory depends on the value 
of the energy// which, in turn, depends on initial conditions. The period of the orbit can 
also be computed as a function of the energy in term of elliptic functions, T = T{H). The 
exact expression is not needed here, but we mention that when the energy tends to zero 
the period tends to infinity, lim//_>o ?"(//) = 00. This periodic movement, see Fig. ([4.^), 
with a well defined period which appears for /i = is characteristic of the KL instability. 

We now interpret the above characteristic of the movement in terms of the non-relaxational 
potential flow plus the orthogonality condition: the relaxational terms in the dynamics make 
the system evolve towards the degenerate minimum of the potential which for fi = occurs 
at oi + 02 + 03 = 1. The residual movement is governed by the non-potential part, see also 



Fig. ( [4.41) . Notice that this residual movement, equations (|4.53D , disappears for 5 = 0, the 
relaxational gradient case. 

Once we have understood the case /i = 0, we now turn to /x > 0. The orthogonality 
condition is no longer satisfied and we can not get such a nice picture of the movement as 
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ai(t),a2(t),a3(t) 

i.G r 




Figure 4.8: Periodic trajectories in the case ^ = 0, 5 = 1.32. The three hnes, solid, dotted 
and dashed, show the time evolution of ai(t), 02 (t) and 03 (t), respectively. Notice that the 
movement is periodic with a well defined period. 



before. However, we can understand what is going on in the following terms: It turns out 
that the energy is no longer a constant of movement after a transient time, but decreases 



exponentially with a characteristics decay time which is of order fi^^ |98]. Consequently, 
and according to the previous analysis, the period of the orbits, which is a function of the 
energy, increases with time. We understand in this way a basic feature of the Busse-Heikes 
model for the KL instability in the case < < 5: the increase of the period between 
successive alternation of the dominating modes, see Fig. ([4.91) . 

This is, indeed, an unwanted feature, since the experimental results do not show such an 
period increase. Busse and Heikes were fully aware of this problem and suggested that noise 
("small amplitude disturbances" as they called them), that is present at all times prevents 
the amplitudes from decaying to arbitrary small levels and a movement which is essentially 
periodic but with a fluctuating period is established.^ We modify, then, Busse-Heikes 
equations by the inclusion of noise terms: 

di = ai(l -ai - 1 +/i(a2 + 03) -(5(a2 -as)) + ^i(t) 

02 = 02(1 - ai - 1 + /x(o2 + 03) - (^("s - «i)) + '^2(*) (4.60) 

03 = 03(1 -ai - 1 +/x(o2 + 03) -(^(fli -«2)) + '^3(*) 



^It is necessary to point out that Tu and Cross ]lO[| have proposed an alternative explanation for the 
stabilization of the period without the necessity of the inclusion of the noise term: they modify the Busse- 
Heike equations by considering a two-dimensional amplitude field and including spatial variation terms. 
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Figure 4.9: Trajectories in the case = 1.20, 5 = 1.32. Notice that the movement is such 
that the time interval between the domination periods of a given mode increases with time. 



where £,i{t) are white-noise processes with correlations: 

mMt')) =2e6{t-t%j (4.61) 

Numerical simulations of these equations for small noise amplitude e shows that the role 
of noise is that of stabilizing the orbits around a mean period. We can understand this in 



the following terms [101]: the inclusion of noise has the effect of injecting "energy" into 
the system. As a consequence, the energy Ti no longer decays to zero but it stabilizes 
around the mean value {H). In this way a periodic movement with a fluctuating period is 



produced, see Fig. ( 4.1C| ). The mean period (T) can be computed from the mean energy 



H by using the same function (T) = T{{H)) that was deduced in the Hamiltonian case. 

We might interpret this in terms of the discussion of Sect. (4.2). For /i = and 5 > 
the function V given by ( |4.43 ) is a Lyapunov potential since it satisfies the orthogonality 



condition. For ^ > this is no longer true but we should expect that for small ^ a 
perturbative solution of the orthogonality condition should give us a potential function 
V^^^ that differs from V it terms that vanish for vanishing fi. In the absence of noise, the 
dynamics leads the system to the minima of the potential V^^^ and the residual movement in 
this attractor is one of increasing period between successive alternation of amplitude modes, 
see Fig. ( [4. 91) . When noise is switched-on, however, fluctuations in the residual motion 
stabilize the mean period to a finite value. The mechanism for this is that fluctuations 
are amplified when the trajectory comes close to one of the unstable fixed points of the 
dynamics and the trajectory is then repelled towards another unstable fixed point. The 
fluctuating period is, hence, sustained by noise. 
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Figure 4.10: Temporal evolution of amplitudes in the case 6 = 1.32, fi = 1.20, e = 10~^. 
In this case, the movement is such that the time interval between dominations of a single 
mode fluctuatues around a mean value. 

5 Noise effects in spatially extended systems 

So far we have considered noise effects in time dependent properties for which spatial degrees 
of freedom are not important. However dynamical noise also has important consequences 



in pattern formation phenomena [95, |102| , 103 1 occurring in spatially extended systems. We 



consider now some of these relevant effects. 

5.1 Symmetry restoring by noise in pattern formation 

Pattern formation is associated with spatial symmetry breaking, but broken symmetries 



might be restored by noise [104]. Equilibrium phase transitions is a well known class of 
problems associated with symmetry breaking. A phase transition takes place when, in 
the thermodynamic limit, thermal fluctuations are not able to mix states with different 
symmetries. This only occurs for a large enough spatial dimensionality, which depends on 
the type of broken symmetry. Noise in phase transitions is discussed in Sect. 6. An example 
of the opposite situation of symmetry restoring by noise is given by the laser instability. 
The description of this instability (Sect. 3) is zero dimensional (that is, described in terms 
of ordinary SDE) and spontaneous emission noise restores the phase symmetry of the lasing 
field after a time of the order of the coherence time: In the long time limit the phase noise 
^(f, implies that there is no preference for any phase <j). Pattern formation usually occurs in 
open systems of low dimensionality, which corresponds to situations somehow intermediate 
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between the two discussed above. Still, in these systems noise may restore the broken 
symmetry implied by the emergence of a pattern, therefore destroying long range order. 
Pattern formation can be sometimes described in terms of ordinary SDE for the amplitudes 
of a few spatial modes, while sometimes a fully continuous description is needed in terms 
of partial SDE. The role of noise in the Kuppers-Lortz instability (Sect. 4) falls within the 
first category, as well as an example of optical pattern formation that we discuss next. A 
prototype of the second category is the stochastic Swift-Hohenberg equation considered later 
in this section. We mention here that the opposite phenomenon of symmetry restoring by 
noise is the one of symmetry breaking by noise in the large system limit. Such phenomenon 
is the one of a noise-induced phase transition considered in Sect. 6. 

The study of spatial and spatiotemporal phenomena in nonlinear optics has emerged 



1 105] as an interesting alternative to more conventional pattern formation studies in fluids. 
For lasers with high Fresnel number (the equivalent of large aspect ratio in hydrodynamics) 
the transverse profile of the laser beam displays pattern formation caused by diffraction 



and nonlinearities. As an example of this, it is known that |10£] when the pump bias level 
of an optically pumped laser is increased above a critical value, a spontaneous breaking of 
cylindrical symmetry in the transverse profile of the laser beam occurs. The emerging state 
displays a spatial pattern associated with phase locking of three appropriate spatial modes. 
The angular position of the pattern fluctuates strongly over time intervals larger than a 
few miliseconds, so that symmetry is restored by a phase diffusion mechanism induced by 
noise ]107|]. We will describe now this phenomenon. Further instabilities of such pattern 



(not considered here) have also been analyzed llOq] . 

The slowly varying envelope E of the electric field inside a ring laser with spherical 
mirrors can be expanded in terms of Gauss-Laguerre cavity modes Ai{p, ip). They describe 
the transverse profile of the beam with p and (p being the polar transverse coordinates. In 
some cavity resonance conditions |107| the electric field can be described in terms of just 
three modes with complex amplitudes /i , /2 , /a : 

3 

E{p,^,t)=Y.mAi{p,^), (5.1) 

i=l 

dth = -fi + 2C[Mifi - ^ Aijkifjfkfn + m, (5.2) 

j,k,l 

where C is the pump parameter, Mi,Aijki given coefficients and model spontaneous 
emission noise as independent Gaussian white noise processes. In this description the 
system becomes effectively zero dimensional with 3 relevant degrees of freedom. Ai is 
independent of (p, while A2, A3 have a cosine and sine profile respectively. In a noise- free 
framework, the laser instability occurs at C = Cth, leading to a cylindrically symmetrical 
pattern with /i 7^ 0, /2 = /s = 0. Increasing further the pump parameter above a second 
threshold C = Ca the modes /2,/3 become unstable and a new pattern emerges. This 
pattern shown in Fig (5.1a) breaks the cylindrical symmetry. The side peaks disappear for 
C < Ca when /2 = /3 = 0. A continuous set of solutions with any angular orientation exists 
and the system spontaneously breaks the symmetry by choosing one of these orientations. 
The dynamical system (5.2) is of the class of relaxational gradient dynamical systems (Sect. 



63 




Figure 5.1: (a) Noise free stationary pattern as obtained from (5.1)-(5.2) and a stationary 
solution of (5.2) for C > Ca- (b) Configuration at time t=150 obtained from the stochastic 
evolution of the configuration in (a) with noise level e = 0.001 

4): 

, fR,i dV R j . , 

where fj^'^ are the real and imaginary parts of /j. Introducing new mode variables gi = 
\9iW^^ : 91 = /i) 52,3 = (/2 =t^/3)/\/2, the angular orientation of the pattern is given by the 
variable r/ = (/?3 — /32)/4. The deterministic stationary solutions of the problem are given 
by the degenerate minima of the potential V for which at I/2I = I/3I, 6 = ±7r/4, ib37r/4, 
with 5 = (/32 + /93)/4 — /3i/2. The discrete values of 5 evidentiate phase locking. 

Symmetry restoring by noise manifests itself in the fact that the potential V , and there- 
fore the stationary solution Pgt (4.36) of the Fokker-Planck equation associated with (5.3) 
is independent of the variable ry. Dynamically, the effect of noise is to restore the cylindrical 
symmetry by driving the system through different values of r/. This effect of noise can be 
visualized by a numerical integration of the stochastic equations for the modal amplitudes 
[|107| ]. Taking as the initial condition a deterministic solution with arbitrary orientation, 
Fig. (5.1a), fluctuations of the modal amplitudes induced by noise lead to a changing pat- 
tern in time with global rotation and distortion. Fig (5.1b). A long time average, or an 
ensemble average over many realizations (as given by Pst)^ results in an averaged cylindri- 
cally symmetric pattern. The dynamics of symmetry restoring is described by the stochastic 
equation for r/: 

dtV = c/i^^ sin45 + i(^e; - ^ei), (5.4) 

where c is a constant proportional to the pump C, /j=|/ip and ^2 3 ^I's real and independent 
Gaussian random processes. The drift part vanishes in the deterministic steady states. For 
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those states, (5.4) describes a diffusive motion of rj (Wiener process) which is the basic 
physical mechanism of symmetry restoring. We also note that the pattern shown in Fig. 
(5.1) has four phase singularities in which E{p,ip) = 0. These singular points also exhibit 



a diffusive motion [107] 



The mechanism of symmetry restoring by phase diffusion also holds in systems which 
require a fully continuous description. If we consider a one dimensional system with a 
stationary pattern forming instability, the relevant pattern forming variable ip can be ex- 
pressed, close to the instability, as [ p^ ] 

V'(x,t) = A(x,t)e^«^^^ + c.c., (5.5) 

where A = Re^^ is the slowly varying amplitude of the most unstable mode qm- Invariance 
under continuous space translations of the original problem [x — > x + xq) implies gauge 
invariance under a change of a constant phase of the amplitude {(p ip + ipQ). Spontaneous 
symmetry breaking associated with pattern formation is here the choice of a particular phase 
of A, which means fixing a global position of the pattern by choosing among a continuous 
family of pattern forming solutions. The phase is the Goldstone mode associated with the 



broken symmetry, which generally satisfies a diffusion equation |10£]: The linear evolution 



of longwavelength modulations of a pattern with dominant wave number q + qM is given by 

dtifqix, t) = D{q)dl^^{x, t) + e(x, t). (5.6) 



The Eckhaus instability | p5| , |110| ] is given by the wavenumber qE for which D{qE) = 0. The 
condition D{q) > identifies the range of wave numbers around qM which give linearly 
stable patterns. Within the range of Eckhaus stable patterns, we have pure stochastic 
diffusion (no drift term) of the A: = Fourier mode component of the phase (pq{x,t). For 
other A;-modes we have damped phase fiuctuations driven by the noise ^. Such stationary 
phase fiuctuations scale as {ip'^{k)) ~ k~'^. The integral of {ip'^{k)) over all wave numbers 
diverges in spatial dimension d < 2. Dimension d = 2 is the critical dimension below which 
fluctuations associated with a broken continuous symmetry destroy long range order, or 
in other words symmetry is restored by noise. In our case this means that fluctuations in 
the local phase of the pattern preclude the existence of a rigidly ordered pattern with well 
defined wavenumber in a large system. We will now be more specific with these ideas in 
the example provided by the Swift-Hohenberg equation. 



5.2 The stochastic Swift-Hohenberg equation: Symmetry restoring, pat- 
tern selection and noise modified Eckhaus instabihty. 

The Swift-Hohenberg equation (SHE) was originally introduced to model the onset of a 
convective instability in simple fiuids |111| . It describes the dynamical evolution of a scalar 
real variable ip{x,t), 

dMx,t) = [/i - (1 + ^2)2] tp{x,t)-tP^{x,t) + V'e^ix,t) = 
5F 

= —^ + V~<{x,t), (5.7) 

where we have explicitly indicated that the SHE has relaxational gradient dynamics de- 
termined by a Lyapunov functional F, and where £,{x, t) models thermal fluctuations as a 
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Figure 5.2: Schematic position of the different wavenumbers mentioned in the text. The 
decay from an Eckhaus unstable wavenumber qi to a final wavenumber Qf is also indicated 

Gaussian process with zero mean and 

t)C{x',t')) = 6{x - x')6{t - t'). (5.8) 
The study of the effect of thermal fluctuations in the onset of pattern formation (^u = 0) was 



one of the original motivations in the study of the SHE [111, 112], but they were commonly 
neglected for some time because of its smallness for simple fluids. However, in several more 
recent experiments, in particular in the electrohydro dynamic instability in nematic liquid 



crystals, these fluctuations have been either directly or indirectly observed [102, 113]. 



In the absence of noise (e = 0) and above threshold (/x > 0), the SHE admits stationary 
periodic solutions tpq{x) characterized by a wavenumber q 

i^qi^) = sm[(2i + l)qx] , q^L < q < QL- (5.9) 

i 

These solutions exist for a band of wavenumbers centered around q = qM = 1 with q±L = 



1 lb These are linearly unstable wave numbers for the no pattern" solution tp = 0, 
being q = qM the most unstable mode. We emphasize that this mode does not coincide 
with the mode qmin associated with the periodic solution which minimizes the potential F. 
Close to threshold one finds [109[ 



qmin ~ 1 - (5-10) 
Stable periodic solutions only exist in a restricted band of wave numbers q-E < q < qE, 



where q±E are the critical wave numbers for the Eckhaus instability described above [114[. 
For small fi, q±E ~ 1 + {q±L — l)/\/3! Fig. (5.2). An important fact about the Eckhaus 
instability in the SHE is that it is a subcritical bifurcation. Defining V{q) = -F[V'g(a;)], 
the Eckhaus boundary is given by the wavenumbers such that ^ = 0, so that qE 
plays the role of a spinodal line in first order phase transitions ]^^ separating unstable 
[qL > q > qE,q-L < q < q-E) from metastable states {q-E < q < qs)- The only globally 
stable state, in the presence of noise, is the one minimizing F. For further reference we also 
note that for the decay of a solution with wavenumber qi outside the Eckhaus stable range, 
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the initially fastest growing mode is determined by a linearized analysis around a periodic 
solution ipq^ix) (Bloch eigenvalue problem). This fastest growing mode is no longer qm = 1 
but qM^Qi), a function of g,. 

We will now address three questions regarding the role of noise in the SHE. First is how 
noise destroys long range order, second the influence of noise in the dynamical process of 
pattern selection from a given initial condition, and finally how noise modifies the Eckhaus 



instability. Answers to these questions 116, 117] are mostly found from extensive 



numerical simulations in which a system of size L, with periodic boundary conditions, 
is discretized in a grid of points Xj, with N ranging from 1024 to 8192 points and a 
discretization step Ax = L/N = 27r/32. For these sizes the number of changes of sign 
(NCS) of a periodic solution with q = qM ranges from 64 to 512, which identifies the 
number of convective rolls in the system. 

Our general arguments above imply that there is no long range order in this one dimen- 
sional system when any amount of noise is present ( a discussion of the d = 2 case is given 
in |ll j j). A quantitative measure of this statement can be given in terms of a stationary 
structure factor defined as 

1 ^ 

hm P(g,t)= |^^(x„f)e-^'?-f . (5.11) 
Alternatively one can look at the normalized correlation function 



{^^J{x + r,t)^Pi x,t)) 



G{r) = lim \^y-^^^'^>^y-^^>l _ (5_^2) 



The numerical results make explicit the absence of long range order in several ways. Con- 
sider first a single realization of the process, that is the stochastic evolution of a configura- 
tion ■0(x,t). The long time dynamics is characterized by mode competition with no single 
dominant mode, so that no periodic solution is established. For small systems (A = 1024) 
this is evidentiated by a power spectrum t) which evolves in time with a hopping max- 
imum wave number and a few modes being permanently excited by noise. The absence 
of a global wave number is also evidentiated by considering a local structure factor at a 
fixed time. Such local structure factor centered in xq is defined replacing in (5.11) 'il^ixi^t) 
by a filtered value around xq: tpxoi^i^'t) — e"^^*"^'"^ il^{xi,t). The local structure factor 
Pxoiq, t), shown in Fig. (5.3), shows a well defined peak at a wave number that is a function 
of xq. Therefore, ip{x,t) is locally periodic, but the local wave number varies from point to 
point. For large systems {N = 8192) the structure factor of a configuration at a fixed time 
is broad, with a typical width of around 40 competing Fourier modes, which gives an idea 
of the range of coherence. In summary, configurations 'ip{x,t) cannot be characterized by a 
single wavenumber. 

Turning now to ensemble averages over several realizations, the stationary correlation 
function shown in Fig. (5.4), is well fitted by 

G(r) = G{0)e~^^^'cos{qMr), (5.13) 

which identifies a correlation length rg. For the typical parameter values of Fig. (5.4) we 
find L/rQ ~ 50 indicating that coherence in the pattern is only maintained in 1/50 of the 
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Figure 5.3: Local power spectrum PxQ{q,t) as defined in the text for a system of N=1024 
points. Parameter values: fi = 0.25, /? = 10 and e = 0.2Aa; 




20 40 60 80 



Figure 5.4: Solid line: Correlation function obtained numerically as an average over 
20 realizations for a system of N=8192 with parameter values /x = 0.25 and e = O.lAx. 
Dashed line: fit to Eq. (5.13) with qm = 1 and ro = 32 
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system size. The Fourier transform of the correlation function, or ensemble averaged power 
spectrum, requires a different fitting for the wings of the spectrum [115] and its width is 
found to grow linearly with noise intensity. We note that despite the absence of long range 
order, and therefore of a wavenumber for the system, still a conservation of the NCS holds 
during the stochastic steady state dynamics. Such conserved NCS can be used to define a 
global wave number or a mean periodicity of the system. 

The question of pattern selection in this problem can be stated as which of the stable 
solutions ipq{x) within the Eckhaus stable band is chosen from a given initial condition. 
Commonly proposed answers to this question are two: a) An argument based on relax- 
ational dynamics on F favor the selection of qmin- This is a potential argument based 
on the role of noise along the dynamical path, b) The argument of fastest linear growth 
favors the selection of qm since this solution is linearly stable, dominates in the early time 
dynamics, and there is conservation of NCS after the initial transient. Two simple com- 
ments are here in order. First is that, from a deterministic point of view, the selected state 
is entirely determined by initial conditions and any linearly stable state is reachable from 
an appropriate initial condition. Secondly is that, from the stochastic point of view, the 
question is ill-posed since we have shown that no single wavenumber characterizes the long 
time configurations. Still, the question makes sense if we ask for the global wave number 
corresponding to a selected NCS. 

From this stochastic point of view one asks for a distribution of selected wavenumbers 
obtained as an ensemble average from a statistically distributed set of typical initial con- 
ditions. A natural situation to consider is the growth of the pattern from random initial 
conditions around ij) = Q. This corresponds again to the case of the decay of an unsta- 
ble state (discussed in Sect. 3), created by changing the control parameter ^ from below 
(/i < 0) to above > 0) the pattern forming instability. In this case the numerical solu- 
tions indicate that noise along path washes out the influence of the initial conflguration. 
The final states are statistically distributed around the NCS corresponding to the global 
wavenumber qm- However, the distribution is significantly narrower than the range of 
Eckhaus stable wavenumbers. Therefore the selected states are not uniformly distributed 
among all possible linearly stable states. The issue of selection of qM versus selection of 
Qmin Cannot be clearly answered numerically because the two wavenumbers are numerically 
very similar. An alternative situation to address this question is to study the decay from a 
solution with wavenumber qi which is Eckhaus unstable. Fig. (5.2). Noise drives the system 
away from the initial condition and there is a transient in which the NCS changes, leading 
to a configuration with a final global wavenumber q^. The time dependence of the NCS for 
a fixed qi and averaged over 20 realizations of the noise is shown in Fig. (5.5). It is seen that 
the evolution gets trapped in a metastable state after the disappearance of 72 rolls. The 
final observed NCS is consistent with the global wavenumber associated with fastest growth 
qAliQi), and it differs from qmin, which minimizes F, in another 36 rolls. Therefore, noise of 
moderate strength along the dynamical path is not able to drive the system, in the explored 
time scale, to the configuration which minimizes the Lyapunov functional of the problem, 
or equivalently to the most probable configuration in terms of the stationary solution of 
the Fokker-Planck equation associated with (5.7). The dynamical mechanism of pattern 
selection seems to be by the mode of fastest growth and this is robust against moderate 
fluctuations. Eventually, and for extremely long time scales, it is a rigorous mathematical 
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Figure 5.5: Time dependence of the averaged NCS in the decay of an Eckhaus unstable 
configuration for a system of N=8192 with fi = 0.56 and e = O.lAx. NCS=512 corresponds 
to g = 1 



statement that noise will induce jumps to configurations with a lower value of F. Of course, 
larger noise intensity would induce earlier jumps, but when the noise becomes larger, the 
conservation of NCS does not hold anymore, and the question of pattern selection becomes 
meaningless. 

The fact that qj is given by the growth of fastest mode qM^Qi) associated with the decay 
from qi breaks down when qi becomes close enough to the Eckhaus stability boundary q±E- 
In this case one enters a fluctuation dominated regime in which qj can be identified neither 
with qhiiqi) nor with qmin- We already mentioned that the Eckhaus boundary is similar to 
a spinodal line in the theory of phase transitions. It is well known in that context that such 
line is only a meanfield concept, so that it is not well defined when fluctuations are present 
A similar results hold here, and the Eckhaus boundary, when fluctuations are present. 



is only defined as an effective stability boundary for a time scale of observation |117]. 

To clarify this point we need first to characterize the range of qi for which we are in a 
fluctuation dominated regime. We can give a criterion which, being similar to the Ginzburg 
criterion of phase transitions [ |119| ], is based on the passage time ideas discussed in Sect. 3. 
The linearized analysis of the dynamics around a configuration ipq. [x) identifies a most un- 
stable Bloch eigenmode qniqi) with eigenvalue or linear growth rate Xm{qi)- The amplitude 
u of this mode satisfies during the early stages of decay a linear SDE 

u{t) = Xmu{t) + V^e(i), (5.14) 

where ^(t) is a white Gaussian random process which results from the projection of ^(x, t) 
onto the most unstable Bloch eigenmode. Equation (5.14), that is the linear regime, is 
valid up to a time T at which the amplitude u becomes large enough, so that nonlinearities 



70 



become important. This time is the MFPT calculated in Sect. 3: 



In(^), (5.15) 



where uq is a reference value for the amplitude u. We should now remember that the 
calculation of T was based on the replacement of a stochastic process h{t), that played 
the role of an effective initial condition, by a stationary random variable h{oo). This 
replacement was valid after a time tm ~ A~^. The time scale given by tm is the fluctuation 
dominated regime in which u{t) stays fluctuating around u = 0. The existence of a linear 
growth regime invoked in our earlier discussion requires a separation of time scales between 
tfn and T. If this is not satisfied we enter directly from a fluctuation dominated regime 
into the nonlinear regime, and the argument of fastest growing mode does not make further 
sense. The separation of time scales tm <C T breaks down for 

Am ~ (5.16) 

Since it is possible to show that Am ~ {qi — Qe)'^, eq. (5.16) determines the range of initial 
wavenumbers qi close to qE for which a linear regime does not exist in the decay process. 
This is the regime dominated by fluctuations for which the definition qe of the Eckhaus 
boundary based on a linear deterministic analysis is not meaningful. A clear example of 
the consequences of this statement is provided by the numerical evidence that, for small 
noise amplitudes, periodic solutions with < qe (within the linear deterministic Eckhaus 
stable band) are not maintained in time and decay to configurations with a global qf closer 
to g = 1. In addition, for initial configurations within the fluctuation dominated regime in 
the Eckhaus unstable range one finds that the result qf = quiqi) does not longer hold. 

Having argued that the Eckhaus boundary is not properly defined in the presence of 
fluctuations, we still can characterize a fluctuation shifted fuzzy boundary qe separating 
unstable from metastable states. Such boundary is defined with respect to a given long 
time scale of observation of the system. Within that scale of observation, it is defined as 
the value of qi for which a relaxation time diverges. Of course this is based on extrapolation 
of data for a range qE < qi in which decay is observed. Numerical results indicate that qf 
becomes a linear function of qi so that qE can be identified by extrapolation to qi = qf. 
We further note that the nonlinear relaxation process in this fluctuation dominated regime 



exhibits a form of dynamical scaling [117|. 



We close this section with a final remark on pattern selection on nonrelaxational systems. 
There are pattern forming models that, as described in Sect. 4, do not follow a relaxational 
dynamics, but still can have a Lyapunov functional. An example recently considered is 



the Greenside-Cross equation in d=2 [123]. Our general discussion on pattern selection 
dynamics makes clear that the configuration minimizing a Lyapunov functional might not 
be the selected configuration in many cases. This has been shown here in the case of a 
purely relaxational dynamics and, a fortiori, will also be true for nonrelaxational dynamics. 
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5.3 Noise amplification in convective instabilities: Noise sustained struc- 
tures 



Generally speaking, we call a state of a system absolutely unstable if a perturbation localized 
around a spatial position grows in time at that spatial position. A state is absolutely stable 
if the perturbation decays in time. An intermediate situation occurs at a convectively 
unstable state: A local perturbation grows in time but traveling in the system, so that 
the perturbation decays in time at fixed space positions but, being convected away from 
the place where it originated, grows in time in the frame of reference that moves with the 



perturbation. The role of noise in a convectively unstable situation is very important |121 
because, as it happened in the decay of an unstable state, there is a mechanism of noise 
amplification: If what grows at the instability is some sort of spatial structure (wave or 
pattern), this structure will not be seen at steady state in the absence of noise because any 
perturbation is convected away and leaves the system. However, if noise is continuously 
present in the system, it will be spatially amplified and a persistent structure sustained 
by noise will be observed. Important questions in this context are the determination of 
conditions for the existence of a convectively unstable regime, the characterization of a noise 
sustained structure and the characterization of the threshold for the transition between a 
noise sustained structure and a deterministically sustained structure. 

As an example of this situation we consider the complex Ginzburg Landau Equation 
(CGLE) which is the amplitude equation for the complex amplitude A of the dominating 
mode at a Hopf bifurcation in a spatially extended system |^5| . Using a different notation 
than in Sect. 4, 

dtA{x, t) - vd^A{x, t) = ijlA{x, t) + {I + ia)dlA{x, t) - {I + il3)\A{x, t)\'^A{x, t) + V^i{x, t) 

(5.17) 

For t; = the state A = changes from absolutely stable to absolutely unstable when the 
control parameter changes from < to ;U > 0. Linearly stable traveling waves exist as 
solutions of the CGLE for > and 1 + a(3 > 0. The convective term vdxA{x, t) becomes 
important when boundary conditions are such that it can not be eliminated by a change of 
frame of reference. To understand the role of this term we linearize (5.17) around A = 0. 
The complex dispersion relation uj for a disturbance of wavenumber K, that thus behaves 
as e'^*"'"^^, becomes: 

oj = n + Kv + {l + ia)K'^ , K = k + iq , (5.18) 

and the growth rate of such a perturbation is given by Reuj{K). Using the method of 
steepest descent, the long-time behavior of the system along a ray defined by fixed x/t, i.e. 
in a frame moving with a velocity vq = x/t, is governed by the saddle point defined by : 

f dLo\ / dio\ 



Re(^^)=v„ , Im{^-)=0. (5.19) 

Since absolute instability occurs when perturbations grow at fixed locations, one has to 
consider the growth rate of modes evolving with zero group velocity, which are defined by: 

/ dio\ , du> 
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These conditions define the following wave number 



q=-ak , k = -^^^^^ . (5.21) 
The real part of uj, which determines the growth rate A of these modes is then: 

Therefore, the uniform reference state {A = 0) is absolutely unstable if A > 0. This 
condition determines a critical line in the parameter space which can be expressed for the 
group velocity v or the control parameter as 



t;^ = 2^/i(l + q2) or = . (5.23) 

Hence, for < // < the uniform reference state is convectively unstable, and wave 
patterns are convected away in the absence of noise. For /i > ^c, wave patterns may grow 
and are sustained by the dynamics, even in the absence of noise. 



This analysis of the convective instability in the CGLE has been used |122] to account 
for the corresponding experimental situation in the case of a Taylor-Couette system with 
through flow, where the transition from convective to absolute instability and noise sus- 
tained structures are observed. Noise sustained structures in the convectively unstable 
regime ^ < fic can be characterized by the stochastic trajectory of yl at a fixed space point. 
It is seen that the statistical properties can be explained in terms of a stochastic phase 
diffusion for the phase of A, while |^| shows small fluctuations. This identifies again a 
structure with no long range order since it is being excited by noise. An equivalent quanti- 
tative measure of the phase diffusion wandering induced by noise amplification is given by 
the width of the frequency power spectrum of A at a fixed point. This spectrum is broad 
in the convectively unstable regime. The threshold at // = //c can be characterized by the 
way in which the width of the spectrum vanishes (divergence of a correlation time) as one 
moves from the convectively unstable regime to the absolute unstable regime. In the latter 
regime the structure persists for vanishing noise and the spectrum is essentially noise free. 



These ideas have been recently extended [124| to the case of two coupled CGLE which 



describe a system undergoing a Hopf bifurcation at a finite wave number. The amplitudes 
A and B of the two emerging counterpropagating waves satisfy 

dtA{x, t) - vd^A{x, t) = nA{x, t) + {I + ia)dlA{x, t) - (1 + i(3)\A{x, t)\'^A{x, t) 

- i^ + i6)\Bix,t)\^Aix,t) + ^Uix,t) (5.24) 
dtB{x,t) +vd^B{x,t) = ^iB{x,t) + {I + ia)dlB{x,t) - {I + if3)\B{x,t)\'^ B{x,t) 

- {J + i6)\A{x,t)fB{x,t) + ^/I^B{x,t). (5.25) 

For 7 > 1 there is strong coupling of the two waves, so that only one of them survives, 
and the situation is the same than for the single CGLE. For 7 < 1 a stability analysis 
analogous to the one discussed for the single CGLE shows that the Traveling Wave solution 
(TW), ^ / 0, i3 = 0, is convectively unstable for fj,c < fJ- < while it becomes absolutely 
unstable for < Therefore, this deterministic analysis would imply that TW solutions 
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Figure 5.6: Space (horizontal axis)-timc (vertical axis) grey scale plot of the moduli of 
the amplitudes A{x,t) (left) and B{x,t) (right) for 7 = 0.8, /j, = 0.408, v = 1, a = 0.02, 
(3 = 0.04, 5 = 0.05. The upper diagrams show the spatial dependence of |^| and \B\ at 
the end of the space-time plots, (a) Deterministic case (e = 0): the disturbances of the 
initial random pattern create initially a standing wave pattern which is replaced, due to 
front propagation, by a symmetric traveling wave pattern, (b) Stochastic case (e = 10~^): 
the spatially distributed noise gives rise, in the bulk of the system, to a noise sustained 
wave structure fluctuating around the deterministic value 
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are observed while they are convectively unstable and Standing Wave solutions (SW), A = 
B ^ 0, emerge when crossing the line 

2 

7c = 1 - - = 1 - TTT^^' (^-26) 

and getting into the region in which the TW solution is absolutely unstable. However, from 
a stochastic point of view one expects that for 7 > 7c (TW convectively unstable) noise will 
be amplified so that the two amplitudes A and B become nonzero and a noise sustained 
SW emerges. Crossing the line 7 = 7c would then lead from a noise sustained SW to a 
deterministically sustained SW. Evidence for the noise sustained SW for 7 > 7c is shown in 
Fig. (5.6), where the results of a numerical simulation with and without noise are shown. 
Initial conditions are random fluctuations around the state A = B = and the boundary 
conditions used in a system of size L are 

A{L,t)=0 , B{0,t)=0 , d^A{0,t)=0 , d^B{L,t) = . (5.27) 

The transition at 7 = 7c can be characterized in terms of the width of the frequency 
power spectrum |j4(3;,C(;)P of A at a fixed space point (inverse correlation time) or by the 
width of the time average of the power spectrum |j4(/c,t)p (inverse correlation length). 
The two power spectra are broad in the region of noise sustained SW. The onset of the 
deterministic SW is identified by a divergence of the correlation time and a correlation 
length becoming of the order of the system size (for the relatively small sizes explored). 
However the transition does not occur at 7 = 7c but at a slightly noise-shifted value |124|. 



This is just another example of instability point shifted and modified by noise as the one 
described above for the Eckhaus instability and the shifting of the phase transition point 
discussed in Sect 6. 

6 Fluctuations, phase transitions and noise— induced transi- 
tions 

In many cases, as we have encountered in previous examples, noise has a disordering ef- 
fect. In this section, however, we consider a class of problem in which noise acts in a non 
trivial and unexpected way inducing a phase transition to an ordered state in spatially 



distributed system |125| , |126| , |127|| . We first review the transitions induced by noise in the 



zero-dimensional case and then turn to the extended systems. 
6.1 Noise— induced transitions 

Noise induced transitions have been known for some time now yi28[ |. Let us start with a 
simple example: consider the following SDE with multiplicative noise (Stratonovich sense) 

± = f{x) + gix)m (6.1) 

^(i) is a white noise process with mean value zero and correlations: 

{^{t)at'))=a''6{t-t') (6.2) 



75 



The parameter a is called the noise intensity and, although it could be absorbed in the 
function g{x), is included here to stress the fact that when o" = the noise effectively 
disappears from the problem. By solving the Fokker-Planck equation, one can check that 
the stationary distribution is given by: 



P.(x)=CexpM (6.3) 
Consequently, the extrema x of the stationary distribution are given by: 

2 

f{x) - '^g{x)g'{x) = (6.4) 

And, quite generally, they will be different from the fixed points of the deterministic dy- 
namics, fj = 0. It is possible, then, that new stationary states appear as a consequence of 
the presence of the noise term. Let us consider a specific case: 

± = -X + \x{l - x^) + {I - x^)i{t) (6.5) 

for < 1 and A < 1. It can be put in the form: 

x = -^ + {l-x^)m (6.6) 

with 

V{x) = 1^x2 + ^x' (6.7) 

In the absence of noise, cr = 0, the deterministic problem has a unique fixed stable point at 
X = 0. For small a the situation is somewhat similar, the stationary probability distribution 
is peaked around the value x = 0. However, by increasing the noise intensity, namely for 
cr^ > 2A, the stationary probability distribution changes shape and becomes bimodal. This 
change from a unimodal to a bimodal distribution obtained by increasing the noise intensity 
has been called a noise-induced transition 0. Of course, this can not be identified as a true 
"phase transition", because during the stochastic evolution the system will jump many 
times from one maximum of the probability to the other, thus not breaking ergodicity, 
which is a clear signature of a phase transition. We show in Fig.( |6.lD how noise is capable 
of restoring the symmetry and, on the average the "order parameter", m = (x) = 0. One 
would think that, in similarity to what happens in other models of statistical mechanics, 
it should be possible to obtain a phase transition by coupling many of these systems. The 
coupling should favor the ordering of neighboring variables to a common value. The system 
would, then, choose between one of the two maxima displaying a macroscopic not vanishing 
order parameter. However, it turns out that, although the noise is effectively capable of 
inducing a phase transition, it does so in a different way. 



A different situation recently studied in [12t] is that of noise induced transition in the hmit of weak 
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t 



Figure 6.1: Typical trajectory coming from a numerical solution of the SDE ( |6.5| ). We see 
that, although there are two preferred values, ±xq, for the variable x (i.e. the stationary 
distribution is bimodal), there are many jumps between these two preferred values. If one 
were to interpret +xo and —xq as two "phases" of the system, we would conclude that noise 
is capable of restoring the symmetry between these two phases. Values of the parameters 
are A = 0.5, ct^ = 2. 



6.2 Noise— induced phase transitions 



Remember that for relaxational systems noise has the effect of introducing fluctuations 
around the minima of the potential. If the potential is bistable and the noise intensity 
becomes large enough, the system is brought out of the minima in a short time scale and 
the order parameter vanishes. One can interpret these results in terms of a shift in the 
transition. Let us be more specific: consider the Ginzburg-Landau model [ 130| , 13l| (model 



A in the taxonomy of |87], also called the (/j^ model). We consider the lattice version: the 
variables Xi follow the relaxational gradient evolution equations: 



dV{x) 
dxi 



(6. 



The index i runs over the 
Ar = 1). The potential is: 



L'^ sites of a d-dimensional regular lattice (of lattice spacing 



N 



v{x) = y: 



1 



1, 



-7:Xi + -X, + -\Vxi 



(6.9) 



i=l 
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We use the simplest form for the lattice gradient Vxj. It is such that its contribution to 
the dynamical equations is the lattice-Laplacian: 



_d_ 

dxi 



2 J2 (^J - ^i) (6-10) 



the sum over the index j runs over the set n{i) of the 2d lattice sites which are nearest 
neighbors of site i. 

As discussed in Sect. 4.2, if the noise term satisfies the fluctuation-dissipation relation: 

{C^{t)Cj{t')) = 2kBTdijd{t - t') (6.11) 

the stationary distribution of (|6.8D is given by: 

Pst{x) = Z-^ exp[-y (x) /ksT] (6.12) 

Z being a normalization constant (the partition function). The parameter h also depends 
on temperature h = h{T) and is such that changes sign at T = Tq, being positive for T < Tq 
and negative for T > Tq. If we neglect thermal fluctuations, the stationary distribution 
becomes a sum of delta functions centered around the minima of the potential V{x). If 
6 < 0, the potential has a single minimum and if 6 > the potential has two minima. In 
these minima all the variables Xj adopt a common value m. The order parameter is: 

Jo 6<0 ,T>ro 

±Vh b>0 ,T<To ^'-'^^ 



This is nothing but a simple example of Landau theory of phase-transitions []132|| . The two 
minima, corresponding to the it signs in the order parameter, are the two phases whose 
symmetry is broken by the time evolution. In the absence of noise, there is no mechanism 
to jump from one minimum to the other. 

If we now go beyond Landau's mean field approximation, we include thermal fluctuations 
by explicitly considering the noise term. Without going into details pH, Il3l|, |l33| , we can 



understand that fluctuations make it possible for the system to jump from one minimum 
to the other, hence restoring the symmetry. For this to happen, the thermal energy ksT 
must be of the order of the energy barrier separating the two minima. Put in another way: 
when thermal fluctuations are taken into account, one needs a deeper potential in order to 
keep the broken symmetry. Or still in other words: when fluctuations are included, there 
is a shift in the location of the phase transition from 6 = to 6 > 0, what implies that 
the new critical temperature is smaller than the Landau mean-field temperature, see Fig. 
(|6.2| ). We conclude that the probability distribution for the order parameter changes from 
bimodal to unimodal when increasing the fluctuations (noise intensity). In fact, for this 
change to happen at temperature greater than zero, one needs the spatial dimension to be 
strictly greater than 1 , d > 1 (lower critical dimension of the f'^ model) . We have witnessed 
in the case of zero-dimensional (one variable) models, how fluctuations can actually go the 
other way and induce a bimodal distribution when increasing the noise intensity. We called 
this a noise-induced transition. In order to have a real noise-induced phase transition, 
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Figure 6.2: Order parameter m for the (p^ model of phase transitions as a function of the 
parameter b. The continuous Une is Landau's solution, eq. ( |6.13| ); it has a critical point at 
6 = 0. The dotted line is the (schematic) exact result (when including thermal fluctuations) 
for spatial dimension d > 1. Notice that there is a shift in the transition when including 
fluctuations towards larger values for 6, smaller temperatures. 



we guess that we need to go a spatially extended system. The role of multiplicative noise 



in spatially extended systems in the context of phase transitions [44, 125, 126, 127 1 and 
pattern-forming instabilities |134|, Q 



136 1 has been thoroughly studied. We consider here 



the d-dimensional version of (|6.1j) with the inclusion of nearest neighbor coupling terms: 



D 



f{xi) + ^ (^i ~ + 9{xi)Ci(t) 



(6.14) 



The parameter D > measures the strength of the coupling between neighboring sites. 
The noise variables ^i(i) are Gaussian distributed of zero mean and correlations given by: 



5{t - 1') 



(6.15) 



We skip now some of the mathematical details. Although it is possible to write down 
the Fokker-Planck equation for the joint probability distribution function P(xi, . . . , xat; i), 
corresponding to this set of stochastic differential equations, it is not possible now, as it 
was in the zero-dimensional case, to solve it in order to find the stationary distribution 
Pst{xi, . . . ,X]\f)- What one does, though, is to obtain an exact equation for the one-site 
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probability Pst{xi): 



_d_ 



d 

-f{xi) + D[xi - {Xn{i)\xi)] + —g{xi)—g{xi) 



Pst{Xi 







(6.16) 



In this equation appears the conditional mean value the mean value of a neighbor 

site of i given that the random variable Xi takes a fixed value. In order to proceed, we make 
now a mean-field type assumption. This approximation, close in spirit to Weiss mean-field 



theory for ferromagnetism [ 13 1 , consists in the assumption that the above conditional mean 
value is equal to the mean value of Xi itself: 



{Xn{i)\Xi) — (Xj) 



(6.17) 



It is possible now to integrate out all the variables except a given one, say xi, from equation 
( 6.16| ) to obtain a closed equation for Pst{xi) whose solution is (we have dropped the 
subindex 1 for the sake of simplicity and because similar equations are satisfied for any Xi): 



Pst{x) 



Z ^ exp 



dy 



f{y)-^9{y)9'{y)-D{y-{x)) 
"^9\y) 



(6.18) 



In this equation, the unknown value (x)is then obtained by the consistency relation: 



dxPst{x) 



(6.19) 



The solutions of equations ( 6.1^ ) and (|6.19|) are to be considered the order parameter, 
m = (x). These equations can be solved numerically, but a first clue about the possible 
solutions can be obtained in the limit D ^ oo by means of a saddle-point integral, which 
yields: 



fi{x)) + -gi{x))g'i{x))=0 



(6.20) 



If by varying the noise intensity, a, the order parameter changes from m = to a non-zero 
value, we interpret this as a phase transition. In the example given below, we will show 
that the main features of a true phase transition hold, namely: the breaking of ergodicity, 
the divergence of correlations, the validity of scaling relations, etc. [|120(| . Finally, note the 
difference in sign of this equation with equation ( |6.4| ) giving the maxima of the probability 
distribution in the case of noise-induced-transitions. This means that, at least in the 
limit of very strong coupling in which equation (6.20) holds, those systems displaying a 
noise-induced transition will not show a noise-induced phase transition and viceversa. 
A simple model is obtained by choosing |125, |126| : 



m 

9{x) 



-x{l + x^f 
1 + 



V{x) = Ul + x 



2^3 



(6.21) 



Let us summarize some properties of this model: 

(i) For the deterministic version (no noise, o" = 0) it is a relaxational gradient system with 
a potential V{x) which displays a single minimum. Therefore, the stationary solution is 
Xi = for all i. 
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Figure 6.3: Numerical solution of equations (|6.18| ) and ( |6.19| ) for the model ( 6.21| ) showing 
the regions of order and disorder (solid line). The dotted line is the result of a higher order 
approximation (not explained in the main text) which goes beyond the approximation ( 6.17] ) 
by including correlations amongst nearest neighbors variables. 



(ii) For the zero-dimensional version, in which there is no coupling between the different 
lattice points, D = 0, each variable becomes independent of the others. By looking at 
the solutions of (|6.4|), one can check that this system does not present a noise-induced 
transition. Therefore, the stationary probability distribution Pst{xi) has a maximum at 
Xj = again for all values of i, and is symmetric around x = 0, hence, m = (xj) = 0. 

(iii) When both coupling and noise terms are present, 7^ 0, cr 7^ 0, the stationary 
distribution has maxima at non-vanishing values of the variable. Due to the symmetry 
X —X, these maxima are at values ±xo, for some xq. In fact, in the infinite coupling 
limit, equation ( |6.2[l| ) yields for the (common) average value: 

m= (x) = ±\/cr^ - 1 (6.22) 

This existence of a non-vanishing order parameter is also maintained at finite values for 
the coupling D, although then the consistency relation ( |6.19| ) must be solved numerically. 
In Fig. ( |6.3D we show the regions in parameter space for which this mean-field approach 
predicts a non-zero value of the order parameter. When the only solution is (x) = we 
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talk about a disordered phase, whereas when non-null solutions (x) ^ appear, we say 
that the system is in an ordered phase. Notice that, according to Fig. (^), an ordered 
phase can appear only for a sufficiently large value of the coupling parameter D and within 
a range of values of the noise intensity a, i.e. the transition is reentrant with respect to the 
noise intensity. 




Figure 6.4: Phase diagram of the system defined by equations ( 6.14| ) and ( 6.21| ) for spatial 



dimension d = 2. The dotted line is the result of the mean field approximation (independent 
of dimension). The symbols are the results of numerical simulations for different system 
sizes (from L = 16, crosses, to L = 128, squares). In this case the order parameter is 
defined as m = (lA^^"*^ J2iLi Xi\)- 



This prediction of the existence of a phase transition in the model of equations ( 6.21| ) 



has been fully confirmed by numerical simulations of the coupled set of Langevin equations 
( 6.14| ) in the two-dimensional, d = 2, case. The numerical solution has used appropriate 



versions of the Milshtein and Heun algorithms, see Sect. (2.4.4). We now describe some 
observed features of the transition as obtained from the computer simulations. 

(i) First of all, there is the breaking of ergodicity. Due to the above mentioned symmetry 
X — > —X and depending on initial conditions and white noise realizations, some realizations 
will go to positive values of the order parameter and some other to negative values. However, 
and this is the key point, the jumps between positive and negative values do not occur very 
often, and the frequency of the jumps strongly decreases with system size N, suggesting 
that in the limit N ^ oo noise is not capable of restoring the symmetry. 

(ii) The second point is the existence of the transition. This is shown in Fig. ( |6.4| ) where 
we plot the order parameter as a function of noise intensity a for a value D = 20 of the 
coupling parameter. According to Fig. (^]^) there should be a non-zero value for the order 
parameter for 1.1 < < 19.9, approximately. Although the reentrant transition occurs at 
a value cr^ ~ 6.7, lower than the predicted in mean field, o"^ ~ 19.9, it is clear the increase 
of the order parameter from m = to a non-zero value and the further decrease towards 
through what appears to be two genuine second-order phase transitions. Of course, calling 
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Figure 6.5: Susceptibility x plotted as a function of the noise intensity . The symbols 
are the results of numerical simulations for different system sizes (same symbols meaning 
as in Fig. ( |6.4D . Notice the sharp increase of the susceptibility as the system size increases, 
as usual in second-order phase transitions. 

this a second-order phase transition is an abuse of language, since there is no free energy 
whose second derivatives show some degree of non-analiticity. Still we keep this name 
because, as we will see in the next paragraphs, many of the characteristics associated to a 
thermodynamic second-order phase transition are also present in this model. 

(iii) In a second-order phase transition, the susceptibility diverges. Here we have defined 
a "susceptibility", X) as the following measure of the fluctuations of the order parameter: 

X = ^r^2 (6-23) 

(the presence of the term is reminiscent of the statistical mechanics relation x = [(w-^) — 
(m)'^]/ [NksT) since in our model the equivalent of temperature T is the noise intensity 
cr^). In Fig. ( |6.5p we can see how the susceptibility develops maxima at the two transitions, 
maxima that increase with increasing system size, suggesting again a true non-analytical 
behavior (divergence) in the thermodynamic limit. It is possible to analyze these results in 
terms of finite size scaling relations, as shown later. 

(iv) Another feature of a second-order phase transition is that of the long-range cor- 
relations. This is shown in Fig. ( |6.6| ) in which we plot the spatial correlation function 
C{n) = ((xj — {xi)){xi^n — i^i+n)))- from the two transition points, the correlation 
function decays very fast, but near both critical points, correlations decay very slowly. The 
same phenomenon has been observed for the correlation times (critical slowing down). 

(v) Finally, we show now that the transition also satisfies the usual scaling relations 



valid in second order transitions. We have focused just in finite-size-scaling relations 137 | 
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n 



Figure 6.6: Correlation function C(n) = ((xj — — for D = 20 and three 

different values of the noise intensity cr^ = 1.65 (solid line), cr^ = 4 (dotted line), = 6.7 
(dashed line). Note that at the two critical points, the correlation function decays very 
slowly. 



(so useful for analysing the results of computer simulations |138, [L39[). Although these 



finite-size scaling relations hold very generally for the model introduced in this section, we 
present results for the order parameter. According to finite size scaling theory, the order 
parameter should behave close enough to the critical points as: 

m(e, L) = L-^/^mieL^/") (6.24) 

where e = 1 — la\ (for thermal phase transitions, this is defined as e = 1 — T/Tc, T being 
the system temperature). The meaning of (3 and v is that of the usual scaling exponents. 
In Fig. (|6.7] ) we plot L^/^m(cr, L) versus 

L^/^(l-crVo-c) to show that curves corresponding 
to different values of L and a fall onto a master curve, so confirming the validity of scaling. 
In this curve we have used the values of the 2-d Ising model exponents, = 1, /? = 1/8. 
Our data, then, suggests that this transition belongs to the 2-d Ising model universality 
class, although more precise data are needed in order to clarify this important point. 

As a conclusion, we have shown in this section how the noise can have a different effect 
that it is usually accepted for Statistical Mechanics systems: it can induce long-range order 
in a spatially extended system. Recent extensions of this work have considered the effect 
of colored noise [ll40(| . This is motivated because one expects that the kind of fluctuations 
leading to multiplicative noise will show a finite correlation time. The role of colored noise 
in the model (with additive noise) has been studied in |141, 142] showing that a non- 
equilibrium phase transition can be induced by varying the correlation time r of the noise. 
For zero-dimensional systems, a reentrant transition has been found in a noise-induced 



transition as a consequence of the color | 143 |. When multiplicative colored noise is included 
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Figure 6.7: Plot of L^/'" m((T, L) versus L^/^{1 — jo^) with a\ = 1.65 for different system 
sizes (same symbols tlian in Fig. (|6.3| )), to show that the finite-size scaHng relation ( 6.24| ) 
is well satisfied. 



in an extended system, several new effects can appear. One of the most counter-intuitive 
ones being that of the disappearance of the order when increasing the coupling between 
nearest neighbor fields |14C| ] . It is hoped that these recent theoretical developments can 
boost experimental work in the search for those predicted phenomena. 
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